From 534d5958810186bb1dee44a4d980e73d89e242d0 Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Fri, 4 Nov 2022 11:39:07 -0700 Subject: [PATCH 01/22] Initial check-in of demo with DMPlex Presently, a DM (of type DMPlex) is created and then information from it is extracted in saved in RDy<*> mesh data structures. --- swe_roe/.gitignore | 1 + swe_roe/ex2.c | 642 +++++++++++++++++++++++++++++++++++++++++++++ swe_roe/makefile | 2 +- 3 files changed, 644 insertions(+), 1 deletion(-) create mode 100644 swe_roe/.gitignore create mode 100644 swe_roe/ex2.c diff --git a/swe_roe/.gitignore b/swe_roe/.gitignore new file mode 100644 index 0000000..cc25a99 --- /dev/null +++ b/swe_roe/.gitignore @@ -0,0 +1 @@ +ex2 diff --git a/swe_roe/ex2.c b/swe_roe/ex2.c new file mode 100644 index 0000000..1ac600d --- /dev/null +++ b/swe_roe/ex2.c @@ -0,0 +1,642 @@ +static char help[] = "Partial 2D dam break problem.\n"; + +#include +#include +#include +#include +#include +#include +#include + +#define RDyAlloc(size, result) PetscCalloc1(size, result) + +typedef struct _n_User *User; + +typedef struct { + PetscReal X[3]; +} RDyCoordinate; + +typedef struct { + PetscReal V[3]; +} RDyVector; + +typedef enum { + CELL_TRI_TYPE = 0, /* tetrahedron cell for a 3D cell */ + CELL_QUAD_TYPE /* hexahedron cell for a 3D cell */ +} RDyCellType; + +typedef struct { + PetscInt *id; /* id of the cell in local numbering */ + PetscInt *global_id; /* global id of the cell in local numbering */ + PetscInt *natural_id; /* natural id of the cell in local numbering */ + + PetscBool *is_local; + + PetscInt *num_vertices; /* number of vertices of the cell */ + PetscInt *num_edges; /* number of edges of the cell */ + PetscInt *num_neighbors; /* number of neigbors of the cell */ + + PetscInt *vertex_offset; /* vertice IDs that form the cell */ + PetscInt *edge_offset; /* vertice IDs that form the cell */ + PetscInt *neighbor_offset; /* vertice IDs that form the cell */ + + PetscInt *vertex_ids; /* vertice IDs that form the cell */ + PetscInt *edge_ids; /* edge IDs that form the cell */ + PetscInt *neighbor_ids; /* neighbor IDs that form the cell */ + + RDyCoordinate *centroid; /* cell centroid */ + + PetscReal *area; /* area of the cell */ + +} RDyCell; + +typedef struct { + PetscInt *id; /* id of the vertex in local numbering */ + PetscInt *global_id; /* global id of the vertex in local numbering */ + + PetscBool *is_local; /* true if the vertex is shared by a local cell */ + + PetscInt *num_cells; /* number of cells sharing the vertex */ + PetscInt *num_edges; /* number of edges sharing the vertex */ + + PetscInt *edge_offset; /* offset for edge IDs that share the vertex */ + PetscInt *cell_offset; /* offset for internal cell IDs that share the vertex */ + + PetscInt *edge_ids; /* edge IDs that share the vertex */ + PetscInt *cell_ids; /* internal cell IDs that share the vertex */ + + RDyCoordinate *coordinate; /* (x,y,z) location of the vertex */ +} RDyVertex; + +typedef struct { + PetscInt *id; /* id of the edge in local numbering */ + PetscInt *global_id; /* global id of the edge in local numbering */ + + PetscBool *is_local; /* true if the edge : (1) */ + /* 1. Is shared by locally owned cells, or */ + /* 2. Is shared by local cell and non-local */ + /* cell such that global ID of local cell */ + /* is smaller than the global ID of */ + /* non-local cell */ + + PetscInt *num_cells; /* number of faces that form the edge */ + PetscInt *vertex_ids; /* ids of vertices that form the edge */ + + PetscInt *cell_offset; /* offset for ids of cell that share the edge */ + PetscInt *cell_ids; /* ids of cells that share the edge */ + + PetscBool *is_internal; /* false if the edge is on the mesh boundary */ + + RDyVector *normal; /* unit normal vector */ + RDyCoordinate *centroid; /* edge centroid */ + + PetscReal *length; /* length of the edge */ + +} RDyEdge; + +typedef struct RDyMesh { + PetscInt dim; + + PetscInt num_cells; + PetscInt num_cells_local; + PetscInt num_edges; + PetscInt num_vertices; + PetscInt num_boundary_faces; + + PetscInt max_vertex_cells, max_vertex_faces; + + RDyCell cells; + RDyVertex vertices; + RDyEdge edges; + + PetscInt *closureSize, **closure, maxClosureSize; + + PetscInt *nG2A; // Mapping of global cells to application/natural cells + +} RDyMesh; + +struct _n_User { + MPI_Comm comm; + char filename[PETSC_MAX_PATH_LEN]; + DM dm; + PetscInt Nt, Nx, Ny; + PetscReal dt, hx, hy; + PetscReal Lx, Ly; + PetscReal hu, hd; + PetscReal tiny_h; + PetscInt dof, rank, comm_size; + Vec B; + PetscBool debug, save, add_building; + PetscInt tstep; + PetscBool interpolate; + + RDyMesh *mesh; +}; + +PetscErrorCode RDyInitialize_IntegerArray_1D(PetscInt *array_1D, PetscInt ndim_1, PetscInt init_value) { + PetscFunctionBegin; + for (PetscInt i = 0; i < ndim_1; i++) { + array_1D[i] = init_value; + } + PetscFunctionReturn(0); +} + +PetscErrorCode RDyAllocate_IntegerArray_1D(PetscInt **array_1D, PetscInt ndim_1) { + PetscFunctionBegin; + PetscCall(RDyAlloc(ndim_1 * sizeof(PetscInt), array_1D)); + PetscCall(RDyInitialize_IntegerArray_1D(*array_1D, ndim_1, -1)); + PetscFunctionReturn(0); +} + +PetscErrorCode RDyAllocate_RDyVector_1D(PetscInt ndim_1, RDyVector **array_1D) { return RDyAlloc(ndim_1 * sizeof(RDyVector), array_1D); } + +PetscErrorCode RDyAllocate_RDyCoordinate_1D(PetscInt ndim_1, RDyCoordinate **array_1D) { return RDyAlloc(ndim_1 * sizeof(RDyCoordinate), array_1D); } + +PetscErrorCode RDyAllocate_RealArray_1D(PetscReal **array_1D, PetscInt ndim_1) { return RDyAlloc(ndim_1 * sizeof(PetscReal), array_1D); } + +PetscBool IsClosureWithinBounds(PetscInt closure, PetscInt start, PetscInt end) { return (closure >= start) && (closure < end); } + +/// @brief Process command line options +/// @param [in] comm A MPI commmunicator +/// @param [inout] user A User data structure that is updated +/// @return 0 on success, or a non-zero error code on failure +static PetscErrorCode ProcessOptions(MPI_Comm comm, User user) { + + PetscFunctionBegin; + + user->comm = comm; + user->Nx = 4; + user->Ny = 5; + user->Lx = user->Nx * 1.0; + user->Ly = user->Ny * 1.0; + user->hx = 1.0; + user->hy = 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; + user->dof = 3; + + MPI_Comm_size(user->comm, &user->comm_size); + MPI_Comm_rank(user->comm, &user->rank); + + PetscErrorCode ierr; + ierr = PetscOptionsBegin(user->comm, NULL, "2D Mesh Options", ""); + PetscCall(ierr); + { + 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("-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(PetscOptionsString("-mesh_filename", "The mesh file", "ex2.c", user->filename, user->filename, PETSC_MAX_PATH_LEN, NULL)); + } + ierr = PetscOptionsEnd(); + + PetscCall(ierr); + assert(user->hu >= 0.); + assert(user->hd >= 0.); + + PetscReal max_time = user->Nt * user->dt; + if (!user->rank) { + PetscPrintf(user->comm, "Max simulation time is %f\n", max_time); + } + + PetscFunctionReturn(0); +} + +/// Creates the PETSc DM as a box or from a file. Add three DOFs and distribute the DM +/// @param [inout] user A User data structure that is modified +/// @return 0 on success, or a non-zero error code on failure +static PetscErrorCode CreateDM(User user) { + + PetscFunctionBegin; + + size_t len; + + PetscStrlen(user->filename, &len); + if (!len) { + PetscInt dim = 2; + PetscInt faces[] = {user->Nx, user->Ny}; + PetscReal lower[] = {0.0, 0.0}; + PetscReal upper[] = {user->Lx, user->Ly}; + + PetscCall(DMPlexCreateBoxMesh(user->comm, dim, PETSC_FALSE, faces, lower, upper, PETSC_NULL, PETSC_TRUE, &user->dm)); + } else { + DMPlexCreateFromFile(user->comm, user->filename, "ex2.c", PETSC_FALSE, &user->dm); + } + + PetscObjectSetName((PetscObject)user->dm, "Mesh"); + PetscCall(DMSetFromOptions(user->dm)); + + // Determine the number of cells, edges, and vertices of the mesh + PetscInt cStart, cEnd; + DMPlexGetHeightStratum(user->dm, 0, &cStart, &cEnd); + + // Create a single section that has 3 DOFs + PetscSection sec; + PetscCall(PetscSectionCreate(user->comm, &sec)); + + // Add the 3 DOFs + PetscInt nfield = 3; + PetscInt num_field_dof[] = {1, 1, 1}; + char field_names[3][20] = {{"Height"}, {"Momentum in x-dir"}, {"Momentum in y-dir"}}; + + nfield = 3; + PetscCall(PetscSectionSetNumFields(sec, nfield)); + PetscInt total_num_dof = 0; + for (PetscInt ifield = 0; ifield < nfield; ifield++) { + PetscCall(PetscSectionSetFieldName(sec, ifield, &field_names[ifield][0])); + PetscCall(PetscSectionSetFieldComponents(sec, ifield, num_field_dof[ifield])); + total_num_dof += num_field_dof[ifield]; + } + + PetscCall(PetscSectionSetChart(sec, cStart, cEnd)); + for (PetscInt c = cStart; c < cEnd; c++) { + for (PetscInt ifield = 0; ifield < nfield; ifield++) { + PetscCall(PetscSectionSetFieldDof(sec, c, ifield, num_field_dof[ifield])); + } + PetscCall(PetscSectionSetDof(sec, c, total_num_dof)); + } + PetscCall(PetscSectionSetUp(sec)); + PetscCall(DMSetLocalSection(user->dm, sec)); + PetscCall(PetscSectionViewFromOptions(sec, NULL, "-layout_view")); + PetscCall(PetscSectionDestroy(&sec)); + PetscCall(DMSetBasicAdjacency(user->dm, PETSC_TRUE, PETSC_TRUE)); + + // Before distributing the DM, set a flag to create mapping from natural-to-local order + PetscCall(DMSetUseNatural(user->dm, PETSC_TRUE)); + + // Distrubte the DM + DM dmDist; + PetscCall(DMPlexDistribute(user->dm, 1, NULL, &dmDist)); + if (dmDist) { + DMDestroy(&user->dm); + user->dm = dmDist; + } + PetscCall(DMViewFromOptions(user->dm, NULL, "-dm_view")); + + PetscFunctionReturn(0); +} + +/// Allocates memory and initialize a RDyCell struct +/// @param [in] num_cells Number of cells +/// @param [out] cells A RDyCell struct that is allocated +/// +/// @return 0 on success, or a non-zero error code on failure +static PetscErrorCode AllocateCells(PetscInt num_cells, RDyCell *cells) { + PetscFunctionBegin; + + PetscInt num_vertices = 4; + PetscInt num_edges = 4; + PetscInt num_neighbors = 4; + + PetscCall(RDyAllocate_IntegerArray_1D(&cells->id, num_cells)); + PetscCall(RDyAllocate_IntegerArray_1D(&cells->global_id, num_cells)); + PetscCall(RDyAllocate_IntegerArray_1D(&cells->natural_id, num_cells)); + + PetscCall(RDyAlloc(num_cells * sizeof(PetscBool), &cells->is_local)); + + PetscCall(RDyAllocate_IntegerArray_1D(&cells->num_vertices, num_cells)); + PetscCall(RDyAllocate_IntegerArray_1D(&cells->num_edges, num_cells)); + PetscCall(RDyAllocate_IntegerArray_1D(&cells->num_neighbors, num_cells)); + + PetscCall(RDyAllocate_IntegerArray_1D(&cells->vertex_offset, num_cells + 1)); + PetscCall(RDyAllocate_IntegerArray_1D(&cells->edge_offset, num_cells + 1)); + PetscCall(RDyAllocate_IntegerArray_1D(&cells->neighbor_offset, num_cells + 1)); + + PetscCall(RDyAllocate_IntegerArray_1D(&cells->vertex_ids, num_cells * num_vertices)); + PetscCall(RDyAllocate_IntegerArray_1D(&cells->edge_ids, num_cells * num_edges)); + PetscCall(RDyAllocate_IntegerArray_1D(&cells->neighbor_ids, num_cells * num_neighbors)); + + PetscCall(RDyAllocate_RDyCoordinate_1D(num_cells, &cells->centroid)); + + PetscCall(RDyAllocate_RealArray_1D(&cells->area, num_cells)); + + for (PetscInt icell = 0; icell < num_cells; icell++) { + cells->id[icell] = icell; + cells->num_vertices[icell] = num_vertices; + cells->num_edges[icell] = num_edges; + cells->num_neighbors[icell] = num_neighbors; + } + + for (PetscInt icell = 0; icell <= num_cells; icell++) { + cells->vertex_offset[icell] = icell * num_vertices; + cells->edge_offset[icell] = icell * num_edges; + cells->neighbor_offset[icell] = icell * num_neighbors; + } + + PetscFunctionReturn(0); +} + +/// Allocates memory and initialize a RDyEdge struct +/// @param [in] num_edges Number of edges +/// @param [out] edges A RDyEdge struct that is allocated +/// +/// @return 0 on success, or a non-zero error code on failure +PetscErrorCode AllocateEdges(PetscInt num_edges, RDyEdge *edges) { + PetscFunctionBegin; + + PetscInt num_cells = 2; + + PetscCall(RDyAllocate_IntegerArray_1D(&edges->id, num_edges)); + PetscCall(RDyAllocate_IntegerArray_1D(&edges->global_id, num_edges)); + PetscCall(RDyAllocate_IntegerArray_1D(&edges->num_cells, num_edges)); + PetscCall(RDyAllocate_IntegerArray_1D(&edges->vertex_ids, num_edges * 2)); + + PetscCall(RDyAlloc(num_edges * sizeof(PetscBool), &edges->is_local)); + PetscCall(RDyAlloc(num_edges * sizeof(PetscBool), &edges->is_internal)); + + PetscCall(RDyAllocate_IntegerArray_1D(&edges->cell_offset, num_edges + 1)); + PetscCall(RDyAllocate_IntegerArray_1D(&edges->cell_ids, num_edges * num_cells)); + + PetscCall(RDyAllocate_RDyCoordinate_1D(num_edges, &edges->centroid)); + PetscCall(RDyAllocate_RDyVector_1D(num_edges, &edges->normal)); + + PetscCall(RDyAllocate_RealArray_1D(&edges->length, num_edges)); + + for (PetscInt iedge = 0; iedge < num_edges; iedge++) { + edges->id[iedge] = iedge; + edges->is_local[iedge] = PETSC_FALSE; + } + for (PetscInt iedge = 0; iedge <= num_edges; iedge++) { + edges->cell_offset[iedge] = iedge * num_cells; + } + + PetscFunctionReturn(0); +} + +/// Allocates memory and initialize a RDyVertex struct. +/// @param [in] num_vertices Number of vertices +/// @param [out] vertices A RDyVertex struct that is allocated +/// +/// @return 0 on success, or a non-zero error code on failure +PetscErrorCode AllocateVertices(PetscInt num_vertices, RDyVertex *vertices) { + PetscFunctionBegin; + + PetscInt ncells_per_vertex = 4; + PetscInt nedges_per_vertex = 4; + + PetscCall(RDyAllocate_IntegerArray_1D(&vertices->id, num_vertices)); + PetscCall(RDyAllocate_IntegerArray_1D(&vertices->global_id, num_vertices)); + PetscCall(RDyAllocate_IntegerArray_1D(&vertices->num_cells, num_vertices)); + PetscCall(RDyAllocate_IntegerArray_1D(&vertices->num_edges, num_vertices)); + + PetscCall(RDyAlloc(num_vertices * sizeof(PetscBool), &vertices->is_local)); + + PetscCall(RDyAllocate_RDyCoordinate_1D(num_vertices, &vertices->coordinate)); + + PetscCall(RDyAllocate_IntegerArray_1D(&vertices->edge_offset, num_vertices + 1)); + PetscCall(RDyAllocate_IntegerArray_1D(&vertices->cell_offset, num_vertices + 1)); + + PetscCall(RDyAllocate_IntegerArray_1D(&vertices->edge_ids, num_vertices * nedges_per_vertex)); + PetscCall(RDyAllocate_IntegerArray_1D(&vertices->cell_ids, num_vertices * ncells_per_vertex)); + + for (PetscInt ivertex = 0; ivertex < num_vertices; ivertex++) { + vertices->id[ivertex] = ivertex; + vertices->is_local[ivertex] = PETSC_FALSE; + vertices->num_cells[ivertex] = 0; + vertices->num_edges[ivertex] = 0; + } + + for (PetscInt ivertex = 0; ivertex <= num_vertices; ivertex++) { + vertices->edge_offset[ivertex] = ivertex * nedges_per_vertex; + vertices->cell_offset[ivertex] = ivertex * ncells_per_vertex; + } + + PetscFunctionReturn(0); +} + +/// @brief Save cell-to-edge, cell-to-vertex, and cell geometric attributes. +/// @param [in] dm A PETSc DM +/// @param [inout] cells A RDyCell structure that is populated +/// @return +static PetscErrorCode PopulateCellsFromDM(DM dm, RDyCell *cells) { + + PetscFunctionBegin; + + PetscInt cStart, cEnd; + PetscInt eStart, eEnd; + DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); + DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd); + + for (PetscInt c = cStart; c < cEnd; c++) { + + PetscInt icell = c - cStart; + PetscInt gref, junkInt; + PetscInt dim=2; + PetscReal centroid[dim], normal[dim]; + PetscCall(DMPlexGetPointGlobal(dm, c, &gref, &junkInt)); + DMPlexComputeCellGeometryFVM(dm, c, &cells->area[icell], ¢roid[0], &normal[0]); + + for (PetscInt idim=0; idimcentroid[icell].X[idim] = centroid[idim]; + } + + PetscInt pSize; + PetscInt *p = NULL; + PetscInt use_cone = PETSC_TRUE; + + cells->num_vertices[icell] = 0; + cells->num_edges[icell] = 0; + if (gref >= 0) { + cells->is_local[icell] = PETSC_TRUE; + } else { + cells->is_local[icell] = PETSC_FALSE; + } + + PetscCall(DMPlexGetTransitiveClosure(dm, c, use_cone, &pSize, &p)); + for (PetscInt i = 2; i < pSize * 2; i += 2) { + if (IsClosureWithinBounds(p[i], eStart, eEnd)) { + PetscInt offset = cells->edge_offset[icell]; + PetscInt index = offset + cells->num_edges[icell]; + cells->edge_ids[index] = p[i]; + cells->num_edges[icell]++; + } else { + PetscInt offset = cells->vertex_offset[icell]; + PetscInt index = offset + cells->num_vertices[icell]; + cells->vertex_ids[index] = p[i]; + cells->num_vertices[icell]++; + } + } + PetscCall(DMPlexRestoreTransitiveClosure(dm, c, use_cone, &pSize, &p)); + } + + PetscFunctionReturn(0); +} + +/// @brief Save edge-to-cell, edge-to-vertex, and geometric attributes. +/// @param [in] dm A PETSc DM +/// @param [inout] edges A RDyVertex structure that is populated +/// @return +static PetscErrorCode PopulateEdgesFromDM(DM dm, RDyEdge *edges) { + + PetscFunctionBegin; + + PetscInt eStart, eEnd; + DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd); + + for (PetscInt e = eStart; e < eEnd; e++) { + + PetscInt iedge = e - eStart; + PetscInt dim = 2; + PetscReal centroid[dim], normal[dim]; + DMPlexComputeCellGeometryFVM(dm, e, &edges->length[iedge], ¢roid[0], &normal[0]); + + for (PetscInt idim=0; idimcentroid[iedge].X[idim] = centroid[idim]; + edges->normal[iedge].V[idim] = normal[idim]; + } + + // edge-to-vertex + PetscInt pSize; + PetscInt *p = NULL; + PetscInt use_cone = PETSC_TRUE; + PetscCall(DMPlexGetTransitiveClosure(dm, e, use_cone, &pSize, &p)); + assert(pSize == 3); + PetscInt index = iedge*2; + edges->vertex_ids[index+0] = p[2]; + edges->vertex_ids[index+1] = p[4]; + PetscCall(DMPlexRestoreTransitiveClosure(dm, e, use_cone, &pSize, &p)); + + // edge-to-cell + edges->num_cells[iedge] = 0; + PetscCall(DMPlexGetTransitiveClosure(dm, e, PETSC_FALSE, &pSize, &p)); + assert(pSize == 2 || pSize == 3); + for (PetscInt i = 2; i < pSize * 2; i += 2) { + PetscInt offset = edges->cell_offset[iedge]; + PetscInt index = offset + edges->num_cells[iedge]; + edges->cell_ids[index] = p[i]; + edges->num_cells[iedge]++; + } + PetscCall(DMPlexRestoreTransitiveClosure(dm, e, PETSC_FALSE, &pSize, &p)); + } + + PetscFunctionReturn(0); +} + +/// @brief Save vertex-to-edge, vertex-to-cell, and geometric attributes +/// (e.g area). +/// @param [in] dm A PETSc DM +/// @param [inout] edges A RDyVertex structure that is populated +/// @return +static PetscErrorCode PopulateVerticesFromDM(DM dm, RDyVertex *vertices) { + + PetscFunctionBegin; + + PetscInt eStart, eEnd; + PetscInt vStart, vEnd; + DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd); + DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd); + + PetscSection coordSection; + Vec coordinates; + DMGetCoordinateSection(dm, &coordSection); + DMGetCoordinatesLocal(dm, &coordinates); + PetscReal *coords; + VecGetArray(coordinates, &coords); + + for (PetscInt v = vStart; v < vEnd; v++) { + + PetscInt ivertex = v - vStart; + PetscInt pSize; + PetscInt *p = NULL; + + PetscCall(DMPlexGetTransitiveClosure(dm, v, PETSC_FALSE, &pSize, &p)); + + PetscInt coordOffset, dim = 2; + PetscSectionGetOffset(coordSection, v, &coordOffset); + for (PetscInt idim=0; idimcoordinate[ivertex].X[idim] = coords[coordOffset+idim]; + } + + vertices->num_edges[ivertex] = 0; + vertices->num_cells[ivertex] = 0; + + for (PetscInt i = 2; i < pSize * 2; i += 2) { + if (IsClosureWithinBounds(p[i], eStart, eEnd)) { + PetscInt offset = vertices->edge_offset[ivertex]; + PetscInt index = offset + vertices->num_edges[ivertex]; + vertices->edge_ids[index] = p[i]; + vertices->num_edges[ivertex]++; + } else { + PetscInt offset = vertices->cell_offset[ivertex]; + PetscInt index = offset + vertices->num_cells[ivertex]; + vertices->cell_ids[index] = p[i]; + vertices->num_cells[ivertex]++; + } + } + + PetscCall(DMPlexRestoreTransitiveClosure(dm, v, PETSC_FALSE, &pSize, &p)); + } + + VecRestoreArray(coordinates, &coords); + + PetscFunctionReturn(0); +} + +/// Creates the RDyMesh structure from PETSc DM +/// @param [inout] user A User data structure that is modified +/// @return 0 on success, or a non-zero error code on failure +static PetscErrorCode CreateMesh(User user) { + + PetscFunctionBegin; + + PetscCall(RDyAlloc(sizeof(RDyMesh), &user->mesh)); + + RDyMesh *mesh_ptr = user->mesh; + + // Determine the number of cells in the mesh + PetscInt cStart, cEnd; + DMPlexGetHeightStratum(user->dm, 0, &cStart, &cEnd); + mesh_ptr->num_cells = cEnd - cStart; + + // Determine the number of edges in the mesh + PetscInt eStart, eEnd; + DMPlexGetDepthStratum(user->dm, 1, &eStart, &eEnd); + mesh_ptr->num_edges = eEnd - eStart; + + // Determine the number of vertices in the mesh + PetscInt vStart, vEnd; + DMPlexGetDepthStratum(user->dm, 0, &vStart, &vEnd); + mesh_ptr->num_vertices = vEnd - vStart; + + // Allocate memory for mesh elements + PetscCall(AllocateCells(mesh_ptr->num_cells, &mesh_ptr->cells)); + PetscCall(AllocateEdges(mesh_ptr->num_edges, &mesh_ptr->edges)); + PetscCall(AllocateVertices(mesh_ptr->num_vertices, &mesh_ptr->vertices)); + + // Populates mesh elements from PETSc DM + PetscCall(PopulateCellsFromDM(user->dm, &mesh_ptr->cells)); + PetscCall(PopulateEdgesFromDM(user->dm, &mesh_ptr->edges)); + PetscCall(PopulateVerticesFromDM(user->dm, &mesh_ptr->vertices)); + + return 0; +} + +int main(int argc, char **argv) { + + PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); + + User user; + PetscCall(PetscNew(&user)); + PetscCall(ProcessOptions(PETSC_COMM_WORLD, user)); + + PetscCall(CreateDM(user)); + + Vec X, R; + PetscCall(DMCreateGlobalVector(user->dm, &X)); // size = dof * number of cells + PetscCall(VecDuplicate(X, &user->B)); + VecViewFromOptions(X, NULL, "-vec_view"); + + PetscCall(CreateMesh(user)); + + PetscCall(PetscFinalize()); + + return 0; +} diff --git a/swe_roe/makefile b/swe_roe/makefile index 3c2e660..b979b56 100644 --- a/swe_roe/makefile +++ b/swe_roe/makefile @@ -4,7 +4,7 @@ FFLAGS = CPPFLAGS = FPPFLAGS = LOCDIR = -EXAMPLESF = ex1.c ex1f.F90 +EXAMPLESF = ex1.c ex1f.F90 ex2.c MANSEC = TS DIRS = toy-problem CLEANFILES = SA-data/* From bd7938669d1e805e57d407d26075a6a88e4c4504 Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Fri, 4 Nov 2022 12:32:40 -0700 Subject: [PATCH 02/22] Adds code to set IC --- swe_roe/ex2.c | 56 ++++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 51 insertions(+), 5 deletions(-) diff --git a/swe_roe/ex2.c b/swe_roe/ex2.c index 1ac600d..9cf2b9f 100644 --- a/swe_roe/ex2.c +++ b/swe_roe/ex2.c @@ -167,8 +167,6 @@ static PetscErrorCode ProcessOptions(MPI_Comm comm, User user) { user->comm = comm; user->Nx = 4; user->Ny = 5; - user->Lx = user->Nx * 1.0; - user->Ly = user->Ny * 1.0; user->hx = 1.0; user->hy = 1.0; user->hu = 10.0; // water depth for the upstream of dam [m] @@ -202,6 +200,9 @@ static PetscErrorCode ProcessOptions(MPI_Comm comm, User user) { assert(user->hu >= 0.); assert(user->hd >= 0.); + user->Lx = user->Nx * 1.0; + user->Ly = user->Ny * 1.0; + PetscReal max_time = user->Nt * user->dt; if (!user->rank) { PetscPrintf(user->comm, "Max simulation time is %f\n", max_time); @@ -416,7 +417,7 @@ PetscErrorCode AllocateVertices(PetscInt num_vertices, RDyVertex *vertices) { /// @param [in] dm A PETSc DM /// @param [inout] cells A RDyCell structure that is populated /// @return -static PetscErrorCode PopulateCellsFromDM(DM dm, RDyCell *cells) { +static PetscErrorCode PopulateCellsFromDM(DM dm, RDyCell *cells, PetscInt *num_cells_local) { PetscFunctionBegin; @@ -425,6 +426,8 @@ static PetscErrorCode PopulateCellsFromDM(DM dm, RDyCell *cells) { DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd); + *num_cells_local = 0; + for (PetscInt c = cStart; c < cEnd; c++) { PetscInt icell = c - cStart; @@ -446,6 +449,7 @@ static PetscErrorCode PopulateCellsFromDM(DM dm, RDyCell *cells) { cells->num_edges[icell] = 0; if (gref >= 0) { cells->is_local[icell] = PETSC_TRUE; + (*num_cells_local)++; } else { cells->is_local[icell] = PETSC_FALSE; } @@ -612,11 +616,39 @@ static PetscErrorCode CreateMesh(User user) { PetscCall(AllocateVertices(mesh_ptr->num_vertices, &mesh_ptr->vertices)); // Populates mesh elements from PETSc DM - PetscCall(PopulateCellsFromDM(user->dm, &mesh_ptr->cells)); + PetscCall(PopulateCellsFromDM(user->dm, &mesh_ptr->cells, &mesh_ptr->num_cells_local)); PetscCall(PopulateEdgesFromDM(user->dm, &mesh_ptr->edges)); PetscCall(PopulateVerticesFromDM(user->dm, &mesh_ptr->vertices)); - return 0; + PetscFunctionReturn(0); +} + + +static PetscErrorCode SetInitialCondition(User user, Vec X) { + + PetscFunctionBegin; + + RDyMesh *mesh = user->mesh; + RDyCell *cells = &mesh->cells; + + PetscCall(VecZeroEntries(X)); + + PetscScalar *x_ptr; + VecGetArray(X, &x_ptr); + + for (PetscInt icell=0; icellnum_cells_local; icell++) { + PetscInt ndof = 3; + PetscInt idx = icell*ndof; + if (cells->centroid[icell].X[1] < 95.0) { + x_ptr[idx] = user->hu; + } else { + x_ptr[idx] = user->hd; + } + } + + VecRestoreArray(X, &x_ptr); + + PetscFunctionReturn(0); } int main(int argc, char **argv) { @@ -627,6 +659,7 @@ int main(int argc, char **argv) { PetscCall(PetscNew(&user)); PetscCall(ProcessOptions(PETSC_COMM_WORLD, user)); + PetscCall(PetscPrintf(PETSC_COMM_WORLD,"1. CreateDM\n")); PetscCall(CreateDM(user)); Vec X, R; @@ -634,8 +667,21 @@ int main(int argc, char **argv) { PetscCall(VecDuplicate(X, &user->B)); VecViewFromOptions(X, NULL, "-vec_view"); + PetscCall(PetscPrintf(PETSC_COMM_WORLD,"2. CreateMesh\n")); PetscCall(CreateMesh(user)); + PetscCall(PetscPrintf(PETSC_COMM_WORLD,"3. SetInitialCondition\n")); + PetscCall(SetInitialCondition(user, X)); + { + char fname[PETSC_MAX_PATH_LEN]; + sprintf(fname, "outputs/ex2_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)); + } + PetscCall(PetscFinalize()); return 0; From 98a634aa9a951f3b3e7b0888bf081e95f3c06239 Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Tue, 8 Nov 2022 08:17:17 -0800 Subject: [PATCH 03/22] Corrects the indices of edges and vertices --- swe_roe/ex2.c | 22 +++++++++++++++------- 1 file changed, 15 insertions(+), 7 deletions(-) diff --git a/swe_roe/ex2.c b/swe_roe/ex2.c index 9cf2b9f..3c83db7 100644 --- a/swe_roe/ex2.c +++ b/swe_roe/ex2.c @@ -423,8 +423,10 @@ static PetscErrorCode PopulateCellsFromDM(DM dm, RDyCell *cells, PetscInt *num_c PetscInt cStart, cEnd; PetscInt eStart, eEnd; + PetscInt vStart, vEnd; DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd); + DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd); *num_cells_local = 0; @@ -459,12 +461,12 @@ static PetscErrorCode PopulateCellsFromDM(DM dm, RDyCell *cells, PetscInt *num_c if (IsClosureWithinBounds(p[i], eStart, eEnd)) { PetscInt offset = cells->edge_offset[icell]; PetscInt index = offset + cells->num_edges[icell]; - cells->edge_ids[index] = p[i]; + cells->edge_ids[index] = p[i] - eStart; cells->num_edges[icell]++; } else { PetscInt offset = cells->vertex_offset[icell]; PetscInt index = offset + cells->num_vertices[icell]; - cells->vertex_ids[index] = p[i]; + cells->vertex_ids[index] = p[i] - vStart; cells->num_vertices[icell]++; } } @@ -482,8 +484,12 @@ static PetscErrorCode PopulateEdgesFromDM(DM dm, RDyEdge *edges) { PetscFunctionBegin; + PetscInt cStart, cEnd; PetscInt eStart, eEnd; + PetscInt vStart, vEnd; + DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd); + DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd); for (PetscInt e = eStart; e < eEnd; e++) { @@ -504,8 +510,8 @@ static PetscErrorCode PopulateEdgesFromDM(DM dm, RDyEdge *edges) { PetscCall(DMPlexGetTransitiveClosure(dm, e, use_cone, &pSize, &p)); assert(pSize == 3); PetscInt index = iedge*2; - edges->vertex_ids[index+0] = p[2]; - edges->vertex_ids[index+1] = p[4]; + edges->vertex_ids[index+0] = p[2] - vStart; + edges->vertex_ids[index+1] = p[4] - vStart; PetscCall(DMPlexRestoreTransitiveClosure(dm, e, use_cone, &pSize, &p)); // edge-to-cell @@ -515,7 +521,7 @@ static PetscErrorCode PopulateEdgesFromDM(DM dm, RDyEdge *edges) { for (PetscInt i = 2; i < pSize * 2; i += 2) { PetscInt offset = edges->cell_offset[iedge]; PetscInt index = offset + edges->num_cells[iedge]; - edges->cell_ids[index] = p[i]; + edges->cell_ids[index] = p[i] - cStart; edges->num_cells[iedge]++; } PetscCall(DMPlexRestoreTransitiveClosure(dm, e, PETSC_FALSE, &pSize, &p)); @@ -533,8 +539,10 @@ static PetscErrorCode PopulateVerticesFromDM(DM dm, RDyVertex *vertices) { PetscFunctionBegin; + PetscInt cStart, cEnd; PetscInt eStart, eEnd; PetscInt vStart, vEnd; + DMPlexGetHeightStratum(dm, 0, &cStart, &cEnd); DMPlexGetDepthStratum(dm, 1, &eStart, &eEnd); DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd); @@ -566,12 +574,12 @@ static PetscErrorCode PopulateVerticesFromDM(DM dm, RDyVertex *vertices) { if (IsClosureWithinBounds(p[i], eStart, eEnd)) { PetscInt offset = vertices->edge_offset[ivertex]; PetscInt index = offset + vertices->num_edges[ivertex]; - vertices->edge_ids[index] = p[i]; + vertices->edge_ids[index] = p[i] - eStart; vertices->num_edges[ivertex]++; } else { PetscInt offset = vertices->cell_offset[ivertex]; PetscInt index = offset + vertices->num_cells[ivertex]; - vertices->cell_ids[index] = p[i]; + vertices->cell_ids[index] = p[i] - cStart; vertices->num_cells[ivertex]++; } } From 746574d8780132bd7c95eaad602b90b2cbebbc19 Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Tue, 8 Nov 2022 09:01:30 -0800 Subject: [PATCH 04/22] Adds TS solver --- swe_roe/ex2.c | 305 +++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 304 insertions(+), 1 deletion(-) diff --git a/swe_roe/ex2.c b/swe_roe/ex2.c index 3c83db7..c75552e 100644 --- a/swe_roe/ex2.c +++ b/swe_roe/ex2.c @@ -126,6 +126,7 @@ struct _n_User { PetscReal tiny_h; PetscInt dof, rank, comm_size; Vec B; + Vec localX; PetscBool debug, save, add_building; PetscInt tstep; PetscBool interpolate; @@ -631,7 +632,10 @@ static PetscErrorCode CreateMesh(User user) { PetscFunctionReturn(0); } - +/// @brief Sets initial condition for [h, hu, hv] +/// @param [in] user A User data structure +/// @param [inout] X Vec for initial condition +/// @return 0 on success, or a non-zero error code on failure static PetscErrorCode SetInitialCondition(User user, Vec X) { PetscFunctionBegin; @@ -659,6 +663,258 @@ static PetscErrorCode SetInitialCondition(User user, Vec X) { PetscFunctionReturn(0); } +PetscErrorCode solver(PetscReal hl, PetscReal hr, PetscReal ul, PetscReal ur, PetscReal vl, PetscReal vr, PetscReal sn, PetscReal cn, + PetscScalar *fij, PetscScalar *amax) { + PetscFunctionBeginUser; + + PetscReal grav = 9.806; + + // Compute Roe averages + PetscReal duml = pow(hl, 0.5); + PetscReal dumr = pow(hr, 0.5); + PetscReal cl = pow(grav * hl, 0.5); + PetscReal cr = pow(grav * hr, 0.5); + PetscReal hhat = duml * dumr; + PetscReal uhat = (duml * ul + dumr * ur) / (duml + dumr); + PetscReal vhat = (duml * vl + dumr * vr) / (duml + dumr); + PetscReal chat = pow(0.5 * grav * (hl + hr), 0.5); + PetscReal uperp = uhat * cn + vhat * sn; + + PetscReal dh = hr - hl; + PetscReal du = ur - ul; + PetscReal dv = vr - vl; + PetscReal dupar = -du * sn + dv * cn; + PetscReal duperp = du * cn + dv * sn; + //printf("cn = %f; sn = %f\n",cn,sn); + + PetscReal dW[3]; + dW[0] = 0.5 * (dh - hhat * duperp / chat); + dW[1] = hhat * dupar; + dW[2] = 0.5 * (dh + hhat * duperp / chat); + + PetscReal uperpl = ul * cn + vl * sn; + PetscReal uperpr = ur * cn + vr * sn; + PetscReal al1 = uperpl - cl; + PetscReal al3 = uperpl + cl; + PetscReal ar1 = uperpr - cr; + PetscReal ar3 = uperpr + cr; + + PetscReal R[3][3]; + R[0][0] = 1.0; + R[0][1] = 0.0; + R[0][2] = 1.0; + R[1][0] = uhat - chat * cn; + R[1][1] = -sn; + R[1][2] = uhat + chat * cn; + R[2][0] = vhat - chat * sn; + R[2][1] = cn; + R[2][2] = vhat + chat * sn; + + PetscReal da1 = fmax(0.0, 2.0 * (ar1 - al1)); + PetscReal da3 = fmax(0.0, 2.0 * (ar3 - al3)); + PetscReal a1 = fabs(uperp - chat); + PetscReal a2 = fabs(uperp); + PetscReal a3 = fabs(uperp + chat); + + // Critical flow fix + if (a1 < da1) { + a1 = 0.5 * (a1 * a1 / da1 + da1); + } + if (a3 < da3) { + a3 = 0.5 * (a3 * a3 / da3 + da3); + } + + // Compute interface flux + PetscReal A[3][3]; + for (PetscInt i = 0; i < 3; i++) { + for (PetscInt j = 0; j < 3; j++) { + A[i][j] = 0.0; + } + } + A[0][0] = a1; + A[1][1] = a2; + A[2][2] = a3; + + PetscReal FL[3], FR[3]; + FL[0] = uperpl * hl; + FL[1] = ul * uperpl * hl + 0.5 * grav * hl * hl * cn; + FL[2] = vl * uperpl * hl + 0.5 * grav * hl * hl * sn; + + FR[0] = uperpr * hr; + FR[1] = ur * uperpr * hr + 0.5 * grav * hr * hr * cn; + FR[2] = vr * uperpr * hr + 0.5 * grav * hr * hr * sn; + + // fij = 0.5*(FL + FR - matmul(R,matmul(A,dW)) + fij[0] = 0.5 * (FL[0] + FR[0] - R[0][0] * A[0][0] * dW[0] - R[0][1] * A[1][1] * dW[1] - R[0][2] * A[2][2] * dW[2]); + fij[1] = 0.5 * (FL[1] + FR[1] - R[1][0] * A[0][0] * dW[0] - R[1][1] * A[1][1] * dW[1] - R[1][2] * A[2][2] * dW[2]); + fij[2] = 0.5 * (FL[2] + FR[2] - R[2][0] * A[0][0] * dW[0] - R[2][1] * A[1][1] * dW[1] - R[2][2] * A[2][2] * dW[2]); + + *amax = chat + fabs(uperp); + + PetscFunctionReturn(0); +} + +static PetscErrorCode GetVelocityFromMomentum(PetscReal tiny_h, PetscReal h, PetscReal hu, PetscReal hv, PetscReal *u, PetscReal *v) { + + PetscFunctionBeginUser; + + if (h < tiny_h) { + *u = 0.0; + *v = 0.0; + } else { + *u = hu/h; + *v = hv/h; + } + + PetscFunctionReturn(0); +} + +PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { + + PetscFunctionBeginUser; + + User user = (User)ptr; + + DM dm = user->dm; + RDyMesh *mesh = user->mesh; + RDyCell *cells = &mesh->cells; + RDyEdge *edges = &mesh->edges; + //RDyVertex *vertices = &mesh->vertices; + + user->tstep = user->tstep + 1; + + PetscCall(DMGlobalToLocalBegin(dm, X, INSERT_VALUES, user->localX)); + PetscCall(DMGlobalToLocalEnd(dm, X, INSERT_VALUES, user->localX)); + PetscCall(VecZeroEntries(F)); + + // Get pointers to vector data + PetscScalar *x_ptr, *f_ptr; + PetscCall(VecGetArray(user->localX, &x_ptr)); + PetscCall(VecGetArray(F, &f_ptr)); + + PetscInt dof = 3; + for (PetscInt iedge=0; iedgenum_edges; iedge++) { + PetscInt cellOffset = edges->cell_offset[iedge]; + PetscInt l = edges->cell_ids[cellOffset ]; + PetscInt r = edges->cell_ids[cellOffset+1]; + + PetscReal sn, cn; + + PetscBool is_edge_vertical = PETSC_TRUE; + if ( PetscAbs(edges->normal[iedge].V[0]) < 1.e-10) { + sn = 1.0; + cn = 0.0; + } else if (PetscAbs(edges->normal[iedge].V[1]) < 1.e-10) { + is_edge_vertical = PETSC_FALSE; + sn = 0.0; + cn = 1.0; + } else { + printf("The code only support quad cells with edges that align with x and y axis\n"); + exit(0); + } + + if (r >= 0 && l >= 0) { + + // Perform computation for an internal edge + + PetscReal hl = x_ptr[l*dof + 0]; + PetscReal hr = x_ptr[r*dof + 0]; + + if (!(hr < user->tiny_h && hl < user->tiny_h)) { + + PetscReal hul = x_ptr[l*dof + 1]; + PetscReal hvl = x_ptr[l*dof + 2]; + PetscReal hur = x_ptr[r*dof + 1]; + PetscReal hvr = x_ptr[r*dof + 2]; + + PetscReal ur, vr, ul, vl; + + PetscCall(GetVelocityFromMomentum(user->tiny_h, hr, hur, hvr, &ur, &vr)); + PetscCall(GetVelocityFromMomentum(user->tiny_h, hl, hul, hvl, &ul, &vl)); + + PetscReal flux[3], amax; + PetscCall(solver(hl, hr, ul, ur, vl, vr, sn, cn, flux, &amax)); + + for (PetscInt idof=0; idofis_local[l]) f_ptr[l*dof + idof] -= flux[idof]; + if (cells->is_local[r]) f_ptr[r*dof + idof] += flux[idof]; + } + } + + } else { + + // Perform computation for a boundary edge + + PetscBool bnd_cell_order_flipped = PETSC_FALSE; + + if (is_edge_vertical) { + if (cells->centroid[l].X[1] > edges->centroid[iedge].X[1]) bnd_cell_order_flipped = PETSC_TRUE; + } else { + if (cells->centroid[l].X[0] > edges->centroid[iedge].X[0]) bnd_cell_order_flipped = PETSC_TRUE; + } + + PetscReal hl = x_ptr[l*dof + 0]; + + if (!(hl < user->tiny_h)) { + + PetscReal hul = x_ptr[l*dof + 1]; + PetscReal hvl = x_ptr[l*dof + 2]; + + PetscReal ul, vl; + PetscCall(GetVelocityFromMomentum(user->tiny_h, hl, hul, hvl, &ul, &vl)); + + PetscReal hr, ur, vr; + hr = hl; + if (is_edge_vertical) { + ur = ul; + vr = -vl; + } else { + ur = -ul; + vr = vl; + } + + if (bnd_cell_order_flipped) { + PetscReal tmp; + tmp = hl; hl = hr; hr = tmp; + tmp = ul; ul = ur; ur = tmp; + tmp = vl; vl = vr; vr = tmp; + } + + PetscReal flux[3], amax; + PetscCall(solver(hl, hr, ul, ur, vl, vr, sn, cn, flux, &amax)); + + for (PetscInt idof=0; idoflocalX, &x_ptr)); + PetscCall(VecRestoreArray(F, &f_ptr)); + + if (user->save) { + char fname[PETSC_MAX_PATH_LEN]; + sprintf(fname, "outputs/ex2_Nx_%d_Ny_%d_dt_%f_%d.dat", user->Nx, user->Ny, user->dt, user->tstep); + PetscViewer viewer; + PetscCall(PetscViewerBinaryOpen(user->comm, fname, FILE_MODE_WRITE, &viewer)); + PetscCall(VecView(X, viewer)); + PetscCall(PetscViewerDestroy(&viewer)); + + sprintf(fname, "outputs/ex2_flux_Nx_%d_Ny_%d_dt_%f_%d.dat", user->Nx, user->Ny, user->dt, user->tstep); + PetscCall(PetscViewerBinaryOpen(user->comm, fname, FILE_MODE_WRITE, &viewer)); + PetscCall(VecView(F, viewer)); + PetscCall(PetscViewerDestroy(&viewer)); + } + + PetscFunctionReturn(0); +} + int main(int argc, char **argv) { PetscCall(PetscInitialize(&argc, &argv, (char *)0, help)); @@ -667,17 +923,31 @@ int main(int argc, char **argv) { PetscCall(PetscNew(&user)); PetscCall(ProcessOptions(PETSC_COMM_WORLD, user)); + /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * + * Create the DM + * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ PetscCall(PetscPrintf(PETSC_COMM_WORLD,"1. CreateDM\n")); PetscCall(CreateDM(user)); + /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * + * Create vectors for solution and residual + * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ Vec X, R; PetscCall(DMCreateGlobalVector(user->dm, &X)); // size = dof * number of cells PetscCall(VecDuplicate(X, &user->B)); + PetscCall(VecDuplicate(X, &R)); VecViewFromOptions(X, NULL, "-vec_view"); + PetscCall(DMCreateLocalVector(user->dm, &user->localX)); // size = dof * number of cells + /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * + * Create the mesh + * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ PetscCall(PetscPrintf(PETSC_COMM_WORLD,"2. CreateMesh\n")); PetscCall(CreateMesh(user)); + /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * + * Initial Condition + * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ PetscCall(PetscPrintf(PETSC_COMM_WORLD,"3. SetInitialCondition\n")); PetscCall(SetInitialCondition(user, X)); { @@ -690,6 +960,39 @@ int main(int argc, char **argv) { PetscCall(PetscViewerDestroy(&viewer)); } + /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * + * Create timestepping solver context + * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + PetscReal max_time = user->Nt * user->dt; + TS ts; + PetscCall(TSCreate(user->comm, &ts)); + PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); + PetscCall(TSSetType(ts, TSEULER)); + PetscCall(TSSetDM(ts, user->dm)); + PetscCall(TSSetRHSFunction(ts, R, RHSFunction, user)); + PetscCall(TSSetMaxTime(ts, max_time)); + PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_STEPOVER)); + PetscCall(TSSetSolution(ts, X)); + PetscCall(TSSetTimeStep(ts, user->dt)); + + PetscCall(TSSetFromOptions(ts)); + PetscCall(TSSolve(ts,X)); + + if (user->save) { + char fname[PETSC_MAX_PATH_LEN]; + sprintf(fname, "outputs/ex2_Nx_%d_Ny_%d_dt_%f_final.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)); + } + + PetscCall(VecDestroy(&X)); + PetscCall(VecDestroy(&R)); + PetscCall(DMDestroy(&user->dm)); + PetscCall(TSDestroy(&ts)); + PetscCall(PetscFinalize()); return 0; From 30c31688153be43071224262618381c33210e8b6 Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Tue, 8 Nov 2022 09:11:53 -0800 Subject: [PATCH 05/22] Adds comments --- swe_roe/ex2.c | 27 +++++++++++++++++++++++++++ 1 file changed, 27 insertions(+) diff --git a/swe_roe/ex2.c b/swe_roe/ex2.c index c75552e..e9b471c 100644 --- a/swe_roe/ex2.c +++ b/swe_roe/ex2.c @@ -663,6 +663,18 @@ static PetscErrorCode SetInitialCondition(User user, Vec X) { PetscFunctionReturn(0); } +/// @brief Computes flux based on Roe solver +/// @param [in] hl Height left of the edge +/// @param [in] hr Height right of the edge +/// @param [in] ul Velocity in x-dir left of the edge +/// @param [in] ur Velocity in x-dir right of the edge +/// @param [in] vl Velocity in y-dir left of the edge +/// @param [in] vr Velocity in y-dir right of the edge +/// @param [in] sn sine of the angle between edge and y-axis +/// @param [in] cn cosine of the angle between edge and y-axis +/// @param [out] fij flux +/// @param [out] amax maximum wave speed +/// @return 0 on success, or a non-zero error code on failure PetscErrorCode solver(PetscReal hl, PetscReal hr, PetscReal ul, PetscReal ur, PetscReal vl, PetscReal vr, PetscReal sn, PetscReal cn, PetscScalar *fij, PetscScalar *amax) { PetscFunctionBeginUser; @@ -754,6 +766,14 @@ PetscErrorCode solver(PetscReal hl, PetscReal hr, PetscReal ul, PetscReal ur, Pe PetscFunctionReturn(0); } +/// @brief Computes velocities in x and y-dir based on momentum in x and y-dir +/// @param tiny_h Threshold value for height +/// @param h Height +/// @param hu Momentum in x-dir +/// @param hv Momentum in y-dir +/// @param u Velocity in x-dir +/// @param v Velocit in y-dir +/// @return 0 on success, or a non-zero error code on failure static PetscErrorCode GetVelocityFromMomentum(PetscReal tiny_h, PetscReal h, PetscReal hu, PetscReal hv, PetscReal *u, PetscReal *v) { PetscFunctionBeginUser; @@ -769,6 +789,13 @@ static PetscErrorCode GetVelocityFromMomentum(PetscReal tiny_h, PetscReal h, Pet PetscFunctionReturn(0); } +/// @brief It is the RHSFunction called by TS +/// @param [in] ts A TS struct +/// @param [in] t Time +/// @param [in] X A global solution Vec +/// @param [inout] F A global flux Vec +/// @param [inout] ptr A user-defined pointer +/// @return 0 on success, or a non-zero error code on failure PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { PetscFunctionBeginUser; From 9bec22c991046b4857d53cafa3d053183fc9217c Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Tue, 8 Nov 2022 10:06:14 -0800 Subject: [PATCH 06/22] Adds code to save cell IDs --- swe_roe/ex2.c | 71 ++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 70 insertions(+), 1 deletion(-) diff --git a/swe_roe/ex2.c b/swe_roe/ex2.c index e9b471c..3dcdf42 100644 --- a/swe_roe/ex2.c +++ b/swe_roe/ex2.c @@ -593,6 +593,74 @@ static PetscErrorCode PopulateVerticesFromDM(DM dm, RDyVertex *vertices) { PetscFunctionReturn(0); } +static PetscErrorCode SaveNaturalCellIDs(DM dm, RDyCell *cells, PetscInt rank) { + + PetscFunctionBegin; + + PetscBool useNatural; + PetscCall(DMGetUseNatural(dm, &useNatural)); + + if (useNatural){ + PetscInt num_fields; + + PetscCall(DMGetNumFields(dm, &num_fields)); + + // Create the natural vector + Vec natural; + PetscCall(DMCreateGlobalVector(dm, &natural)); + PetscInt natural_size = 0, cum_natural_size = 0; + PetscCall(VecGetLocalSize(natural, &natural_size)); + PetscCall(MPI_Scan(&natural_size, &cum_natural_size, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD)); + + // Add entries in the natural vector + PetscScalar *entries; + printf("rank = %d; num_fields = %d; natural_size = %d %d; offset = %d\n",rank, num_fields,natural_size,cum_natural_size, cum_natural_size/num_fields - natural_size/num_fields); + PetscCall(VecGetArray(natural, &entries)); + for (PetscInt i = 0; i < natural_size; ++i) { + if (i % num_fields == 0) { + entries[i] = i/num_fields + cum_natural_size/num_fields - natural_size/num_fields; + } + else { + entries[i] = -1 - rank; + } + } + PetscCall(VecRestoreArray(natural, &entries)); + VecView(natural,PETSC_VIEWER_STDOUT_WORLD); + + // Map natural IDs in global order + Vec global; + PetscCall(DMCreateGlobalVector(dm, &global)); + PetscCall(DMPlexNaturalToGlobalBegin(dm, natural, global)); + PetscCall(DMPlexNaturalToGlobalEnd(dm, natural, global)); + VecView(global,PETSC_VIEWER_STDOUT_WORLD); + + // Map natural IDs in local order + Vec local; + PetscCall(DMCreateLocalVector(dm, &local)); + PetscCall(DMGlobalToLocalBegin(dm, global, INSERT_VALUES, local)); + PetscCall(DMGlobalToLocalEnd(dm, global, INSERT_VALUES, local)); + VecView(local,PETSC_VIEWER_STDOUT_WORLD); + + // Save natural IDs + PetscInt local_size; + PetscCall(VecGetLocalSize(local, &local_size)); + PetscCall(VecGetArray(local, &entries)); + printf("rank = %d; local_size = %d; \n",rank,local_size); + for (PetscInt i = 0; i < local_size/num_fields; ++i) { + cells->natural_id[i] = entries[i*num_fields]; + } + PetscCall(VecRestoreArray(local, &entries)); + + // Cleanup + PetscCall(VecDestroy(&natural)); + PetscCall(VecDestroy(&global)); + PetscCall(VecDestroy(&local)); + + } + + PetscFunctionReturn(0); +} + /// Creates the RDyMesh structure from PETSc DM /// @param [inout] user A User data structure that is modified /// @return 0 on success, or a non-zero error code on failure @@ -628,6 +696,7 @@ static PetscErrorCode CreateMesh(User user) { PetscCall(PopulateCellsFromDM(user->dm, &mesh_ptr->cells, &mesh_ptr->num_cells_local)); PetscCall(PopulateEdgesFromDM(user->dm, &mesh_ptr->edges)); PetscCall(PopulateVerticesFromDM(user->dm, &mesh_ptr->vertices)); + PetscCall(SaveNaturalCellIDs(user->dm, &mesh_ptr->cells, user->rank)); PetscFunctionReturn(0); } @@ -868,7 +937,7 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { } } - } else { + } else if (cells->is_local[l]) { // Perform computation for a boundary edge From 9d9b79404e7ce3f4f5b24a2472121979a5c5db0f Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Tue, 8 Nov 2022 10:23:26 -0800 Subject: [PATCH 07/22] Fixes formatting --- swe_roe/ex2.c | 195 +++++++++++++++++++++++--------------------------- 1 file changed, 91 insertions(+), 104 deletions(-) diff --git a/swe_roe/ex2.c b/swe_roe/ex2.c index 3dcdf42..58884d6 100644 --- a/swe_roe/ex2.c +++ b/swe_roe/ex2.c @@ -57,12 +57,12 @@ typedef struct { PetscBool *is_local; /* true if the vertex is shared by a local cell */ PetscInt *num_cells; /* number of cells sharing the vertex */ - PetscInt *num_edges; /* number of edges sharing the vertex */ + PetscInt *num_edges; /* number of edges sharing the vertex */ - PetscInt *edge_offset; /* offset for edge IDs that share the vertex */ + PetscInt *edge_offset; /* offset for edge IDs that share the vertex */ PetscInt *cell_offset; /* offset for internal cell IDs that share the vertex */ - PetscInt *edge_ids; /* edge IDs that share the vertex */ + PetscInt *edge_ids; /* edge IDs that share the vertex */ PetscInt *cell_ids; /* internal cell IDs that share the vertex */ RDyCoordinate *coordinate; /* (x,y,z) location of the vertex */ @@ -162,7 +162,6 @@ PetscBool IsClosureWithinBounds(PetscInt closure, PetscInt start, PetscInt end) /// @param [inout] user A User data structure that is updated /// @return 0 on success, or a non-zero error code on failure static PetscErrorCode ProcessOptions(MPI_Comm comm, User user) { - PetscFunctionBegin; user->comm = comm; @@ -216,7 +215,6 @@ static PetscErrorCode ProcessOptions(MPI_Comm comm, User user) { /// @param [inout] user A User data structure that is modified /// @return 0 on success, or a non-zero error code on failure static PetscErrorCode CreateDM(User user) { - PetscFunctionBegin; size_t len; @@ -416,10 +414,9 @@ PetscErrorCode AllocateVertices(PetscInt num_vertices, RDyVertex *vertices) { /// @brief Save cell-to-edge, cell-to-vertex, and cell geometric attributes. /// @param [in] dm A PETSc DM -/// @param [inout] cells A RDyCell structure that is populated -/// @return +/// @param [inout] cells A RDyCell structure that is populated +/// @return static PetscErrorCode PopulateCellsFromDM(DM dm, RDyCell *cells, PetscInt *num_cells_local) { - PetscFunctionBegin; PetscInt cStart, cEnd; @@ -432,15 +429,14 @@ static PetscErrorCode PopulateCellsFromDM(DM dm, RDyCell *cells, PetscInt *num_c *num_cells_local = 0; for (PetscInt c = cStart; c < cEnd; c++) { - - PetscInt icell = c - cStart; + PetscInt icell = c - cStart; PetscInt gref, junkInt; - PetscInt dim=2; + PetscInt dim = 2; PetscReal centroid[dim], normal[dim]; PetscCall(DMPlexGetPointGlobal(dm, c, &gref, &junkInt)); DMPlexComputeCellGeometryFVM(dm, c, &cells->area[icell], ¢roid[0], &normal[0]); - for (PetscInt idim=0; idimcentroid[icell].X[idim] = centroid[idim]; } @@ -449,7 +445,7 @@ static PetscErrorCode PopulateCellsFromDM(DM dm, RDyCell *cells, PetscInt *num_c PetscInt use_cone = PETSC_TRUE; cells->num_vertices[icell] = 0; - cells->num_edges[icell] = 0; + cells->num_edges[icell] = 0; if (gref >= 0) { cells->is_local[icell] = PETSC_TRUE; (*num_cells_local)++; @@ -460,13 +456,13 @@ static PetscErrorCode PopulateCellsFromDM(DM dm, RDyCell *cells, PetscInt *num_c PetscCall(DMPlexGetTransitiveClosure(dm, c, use_cone, &pSize, &p)); for (PetscInt i = 2; i < pSize * 2; i += 2) { if (IsClosureWithinBounds(p[i], eStart, eEnd)) { - PetscInt offset = cells->edge_offset[icell]; - PetscInt index = offset + cells->num_edges[icell]; + PetscInt offset = cells->edge_offset[icell]; + PetscInt index = offset + cells->num_edges[icell]; cells->edge_ids[index] = p[i] - eStart; cells->num_edges[icell]++; } else { - PetscInt offset = cells->vertex_offset[icell]; - PetscInt index = offset + cells->num_vertices[icell]; + PetscInt offset = cells->vertex_offset[icell]; + PetscInt index = offset + cells->num_vertices[icell]; cells->vertex_ids[index] = p[i] - vStart; cells->num_vertices[icell]++; } @@ -479,10 +475,9 @@ static PetscErrorCode PopulateCellsFromDM(DM dm, RDyCell *cells, PetscInt *num_c /// @brief Save edge-to-cell, edge-to-vertex, and geometric attributes. /// @param [in] dm A PETSc DM -/// @param [inout] edges A RDyVertex structure that is populated -/// @return +/// @param [inout] edges A RDyVertex structure that is populated +/// @return static PetscErrorCode PopulateEdgesFromDM(DM dm, RDyEdge *edges) { - PetscFunctionBegin; PetscInt cStart, cEnd; @@ -493,15 +488,14 @@ static PetscErrorCode PopulateEdgesFromDM(DM dm, RDyEdge *edges) { DMPlexGetDepthStratum(dm, 0, &vStart, &vEnd); for (PetscInt e = eStart; e < eEnd; e++) { - - PetscInt iedge = e - eStart; - PetscInt dim = 2; + PetscInt iedge = e - eStart; + PetscInt dim = 2; PetscReal centroid[dim], normal[dim]; DMPlexComputeCellGeometryFVM(dm, e, &edges->length[iedge], ¢roid[0], &normal[0]); - for (PetscInt idim=0; idimcentroid[iedge].X[idim] = centroid[idim]; - edges->normal[iedge].V[idim] = normal[idim]; + edges->normal[iedge].V[idim] = normal[idim]; } // edge-to-vertex @@ -510,9 +504,9 @@ static PetscErrorCode PopulateEdgesFromDM(DM dm, RDyEdge *edges) { PetscInt use_cone = PETSC_TRUE; PetscCall(DMPlexGetTransitiveClosure(dm, e, use_cone, &pSize, &p)); assert(pSize == 3); - PetscInt index = iedge*2; - edges->vertex_ids[index+0] = p[2] - vStart; - edges->vertex_ids[index+1] = p[4] - vStart; + PetscInt index = iedge * 2; + edges->vertex_ids[index + 0] = p[2] - vStart; + edges->vertex_ids[index + 1] = p[4] - vStart; PetscCall(DMPlexRestoreTransitiveClosure(dm, e, use_cone, &pSize, &p)); // edge-to-cell @@ -520,8 +514,8 @@ static PetscErrorCode PopulateEdgesFromDM(DM dm, RDyEdge *edges) { PetscCall(DMPlexGetTransitiveClosure(dm, e, PETSC_FALSE, &pSize, &p)); assert(pSize == 2 || pSize == 3); for (PetscInt i = 2; i < pSize * 2; i += 2) { - PetscInt offset = edges->cell_offset[iedge]; - PetscInt index = offset + edges->num_cells[iedge]; + PetscInt offset = edges->cell_offset[iedge]; + PetscInt index = offset + edges->num_cells[iedge]; edges->cell_ids[index] = p[i] - cStart; edges->num_cells[iedge]++; } @@ -531,13 +525,12 @@ static PetscErrorCode PopulateEdgesFromDM(DM dm, RDyEdge *edges) { PetscFunctionReturn(0); } -/// @brief Save vertex-to-edge, vertex-to-cell, and geometric attributes +/// @brief Save vertex-to-edge, vertex-to-cell, and geometric attributes /// (e.g area). /// @param [in] dm A PETSc DM -/// @param [inout] edges A RDyVertex structure that is populated -/// @return +/// @param [inout] edges A RDyVertex structure that is populated +/// @return static PetscErrorCode PopulateVerticesFromDM(DM dm, RDyVertex *vertices) { - PetscFunctionBegin; PetscInt cStart, cEnd; @@ -555,8 +548,7 @@ static PetscErrorCode PopulateVerticesFromDM(DM dm, RDyVertex *vertices) { VecGetArray(coordinates, &coords); for (PetscInt v = vStart; v < vEnd; v++) { - - PetscInt ivertex = v - vStart; + PetscInt ivertex = v - vStart; PetscInt pSize; PetscInt *p = NULL; @@ -564,8 +556,8 @@ static PetscErrorCode PopulateVerticesFromDM(DM dm, RDyVertex *vertices) { PetscInt coordOffset, dim = 2; PetscSectionGetOffset(coordSection, v, &coordOffset); - for (PetscInt idim=0; idimcoordinate[ivertex].X[idim] = coords[coordOffset+idim]; + for (PetscInt idim = 0; idim < dim; idim++) { + vertices->coordinate[ivertex].X[idim] = coords[coordOffset + idim]; } vertices->num_edges[ivertex] = 0; @@ -573,13 +565,13 @@ static PetscErrorCode PopulateVerticesFromDM(DM dm, RDyVertex *vertices) { for (PetscInt i = 2; i < pSize * 2; i += 2) { if (IsClosureWithinBounds(p[i], eStart, eEnd)) { - PetscInt offset = vertices->edge_offset[ivertex]; - PetscInt index = offset + vertices->num_edges[ivertex]; + PetscInt offset = vertices->edge_offset[ivertex]; + PetscInt index = offset + vertices->num_edges[ivertex]; vertices->edge_ids[index] = p[i] - eStart; vertices->num_edges[ivertex]++; } else { - PetscInt offset = vertices->cell_offset[ivertex]; - PetscInt index = offset + vertices->num_cells[ivertex]; + PetscInt offset = vertices->cell_offset[ivertex]; + PetscInt index = offset + vertices->num_cells[ivertex]; vertices->cell_ids[index] = p[i] - cStart; vertices->num_cells[ivertex]++; } @@ -594,13 +586,12 @@ static PetscErrorCode PopulateVerticesFromDM(DM dm, RDyVertex *vertices) { } static PetscErrorCode SaveNaturalCellIDs(DM dm, RDyCell *cells, PetscInt rank) { - PetscFunctionBegin; PetscBool useNatural; PetscCall(DMGetUseNatural(dm, &useNatural)); - if (useNatural){ + if (useNatural) { PetscInt num_fields; PetscCall(DMGetNumFields(dm, &num_fields)); @@ -614,40 +605,40 @@ static PetscErrorCode SaveNaturalCellIDs(DM dm, RDyCell *cells, PetscInt rank) { // Add entries in the natural vector PetscScalar *entries; - printf("rank = %d; num_fields = %d; natural_size = %d %d; offset = %d\n",rank, num_fields,natural_size,cum_natural_size, cum_natural_size/num_fields - natural_size/num_fields); + printf("rank = %d; num_fields = %d; natural_size = %d %d; offset = %d\n", rank, num_fields, natural_size, cum_natural_size, + cum_natural_size / num_fields - natural_size / num_fields); PetscCall(VecGetArray(natural, &entries)); for (PetscInt i = 0; i < natural_size; ++i) { if (i % num_fields == 0) { - entries[i] = i/num_fields + cum_natural_size/num_fields - natural_size/num_fields; - } - else { + entries[i] = i / num_fields + cum_natural_size / num_fields - natural_size / num_fields; + } else { entries[i] = -1 - rank; } } PetscCall(VecRestoreArray(natural, &entries)); - VecView(natural,PETSC_VIEWER_STDOUT_WORLD); + VecView(natural, PETSC_VIEWER_STDOUT_WORLD); // Map natural IDs in global order Vec global; PetscCall(DMCreateGlobalVector(dm, &global)); PetscCall(DMPlexNaturalToGlobalBegin(dm, natural, global)); PetscCall(DMPlexNaturalToGlobalEnd(dm, natural, global)); - VecView(global,PETSC_VIEWER_STDOUT_WORLD); + VecView(global, PETSC_VIEWER_STDOUT_WORLD); // Map natural IDs in local order Vec local; PetscCall(DMCreateLocalVector(dm, &local)); PetscCall(DMGlobalToLocalBegin(dm, global, INSERT_VALUES, local)); PetscCall(DMGlobalToLocalEnd(dm, global, INSERT_VALUES, local)); - VecView(local,PETSC_VIEWER_STDOUT_WORLD); + VecView(local, PETSC_VIEWER_STDOUT_WORLD); // Save natural IDs PetscInt local_size; PetscCall(VecGetLocalSize(local, &local_size)); PetscCall(VecGetArray(local, &entries)); - printf("rank = %d; local_size = %d; \n",rank,local_size); - for (PetscInt i = 0; i < local_size/num_fields; ++i) { - cells->natural_id[i] = entries[i*num_fields]; + printf("rank = %d; local_size = %d; \n", rank, local_size); + for (PetscInt i = 0; i < local_size / num_fields; ++i) { + cells->natural_id[i] = entries[i * num_fields]; } PetscCall(VecRestoreArray(local, &entries)); @@ -655,7 +646,6 @@ static PetscErrorCode SaveNaturalCellIDs(DM dm, RDyCell *cells, PetscInt rank) { PetscCall(VecDestroy(&natural)); PetscCall(VecDestroy(&global)); PetscCall(VecDestroy(&local)); - } PetscFunctionReturn(0); @@ -665,7 +655,6 @@ static PetscErrorCode SaveNaturalCellIDs(DM dm, RDyCell *cells, PetscInt rank) { /// @param [inout] user A User data structure that is modified /// @return 0 on success, or a non-zero error code on failure static PetscErrorCode CreateMesh(User user) { - PetscFunctionBegin; PetscCall(RDyAlloc(sizeof(RDyMesh), &user->mesh)); @@ -706,10 +695,9 @@ static PetscErrorCode CreateMesh(User user) { /// @param [inout] X Vec for initial condition /// @return 0 on success, or a non-zero error code on failure static PetscErrorCode SetInitialCondition(User user, Vec X) { - PetscFunctionBegin; - RDyMesh *mesh = user->mesh; + RDyMesh *mesh = user->mesh; RDyCell *cells = &mesh->cells; PetscCall(VecZeroEntries(X)); @@ -717,9 +705,9 @@ static PetscErrorCode SetInitialCondition(User user, Vec X) { PetscScalar *x_ptr; VecGetArray(X, &x_ptr); - for (PetscInt icell=0; icellnum_cells_local; icell++) { + for (PetscInt icell = 0; icell < mesh->num_cells_local; icell++) { PetscInt ndof = 3; - PetscInt idx = icell*ndof; + PetscInt idx = icell * ndof; if (cells->centroid[icell].X[1] < 95.0) { x_ptr[idx] = user->hu; } else { @@ -766,7 +754,7 @@ PetscErrorCode solver(PetscReal hl, PetscReal hr, PetscReal ul, PetscReal ur, Pe PetscReal dv = vr - vl; PetscReal dupar = -du * sn + dv * cn; PetscReal duperp = du * cn + dv * sn; - //printf("cn = %f; sn = %f\n",cn,sn); + // printf("cn = %f; sn = %f\n",cn,sn); PetscReal dW[3]; dW[0] = 0.5 * (dh - hhat * duperp / chat); @@ -844,15 +832,14 @@ PetscErrorCode solver(PetscReal hl, PetscReal hr, PetscReal ul, PetscReal ur, Pe /// @param v Velocit in y-dir /// @return 0 on success, or a non-zero error code on failure static PetscErrorCode GetVelocityFromMomentum(PetscReal tiny_h, PetscReal h, PetscReal hu, PetscReal hv, PetscReal *u, PetscReal *v) { - PetscFunctionBeginUser; if (h < tiny_h) { *u = 0.0; *v = 0.0; } else { - *u = hu/h; - *v = hv/h; + *u = hu / h; + *v = hv / h; } PetscFunctionReturn(0); @@ -866,16 +853,15 @@ static PetscErrorCode GetVelocityFromMomentum(PetscReal tiny_h, PetscReal h, Pet /// @param [inout] ptr A user-defined pointer /// @return 0 on success, or a non-zero error code on failure PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { - PetscFunctionBeginUser; User user = (User)ptr; - DM dm = user->dm; - RDyMesh *mesh = user->mesh; + DM dm = user->dm; + RDyMesh *mesh = user->mesh; RDyCell *cells = &mesh->cells; RDyEdge *edges = &mesh->edges; - //RDyVertex *vertices = &mesh->vertices; + // RDyVertex *vertices = &mesh->vertices; user->tstep = user->tstep + 1; @@ -889,39 +875,37 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { PetscCall(VecGetArray(F, &f_ptr)); PetscInt dof = 3; - for (PetscInt iedge=0; iedgenum_edges; iedge++) { + for (PetscInt iedge = 0; iedge < mesh->num_edges; iedge++) { PetscInt cellOffset = edges->cell_offset[iedge]; - PetscInt l = edges->cell_ids[cellOffset ]; - PetscInt r = edges->cell_ids[cellOffset+1]; + PetscInt l = edges->cell_ids[cellOffset]; + PetscInt r = edges->cell_ids[cellOffset + 1]; PetscReal sn, cn; PetscBool is_edge_vertical = PETSC_TRUE; - if ( PetscAbs(edges->normal[iedge].V[0]) < 1.e-10) { + if (PetscAbs(edges->normal[iedge].V[0]) < 1.e-10) { sn = 1.0; cn = 0.0; } else if (PetscAbs(edges->normal[iedge].V[1]) < 1.e-10) { is_edge_vertical = PETSC_FALSE; - sn = 0.0; - cn = 1.0; + sn = 0.0; + cn = 1.0; } else { printf("The code only support quad cells with edges that align with x and y axis\n"); exit(0); } if (r >= 0 && l >= 0) { - // Perform computation for an internal edge - PetscReal hl = x_ptr[l*dof + 0]; - PetscReal hr = x_ptr[r*dof + 0]; + PetscReal hl = x_ptr[l * dof + 0]; + PetscReal hr = x_ptr[r * dof + 0]; if (!(hr < user->tiny_h && hl < user->tiny_h)) { - - PetscReal hul = x_ptr[l*dof + 1]; - PetscReal hvl = x_ptr[l*dof + 2]; - PetscReal hur = x_ptr[r*dof + 1]; - PetscReal hvr = x_ptr[r*dof + 2]; + PetscReal hul = x_ptr[l * dof + 1]; + PetscReal hvl = x_ptr[l * dof + 2]; + PetscReal hur = x_ptr[r * dof + 1]; + PetscReal hvr = x_ptr[r * dof + 2]; PetscReal ur, vr, ul, vl; @@ -931,14 +915,13 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { PetscReal flux[3], amax; PetscCall(solver(hl, hr, ul, ur, vl, vr, sn, cn, flux, &amax)); - for (PetscInt idof=0; idofis_local[l]) f_ptr[l*dof + idof] -= flux[idof]; - if (cells->is_local[r]) f_ptr[r*dof + idof] += flux[idof]; + for (PetscInt idof = 0; idof < dof; idof++) { + if (cells->is_local[l]) f_ptr[l * dof + idof] -= flux[idof]; + if (cells->is_local[r]) f_ptr[r * dof + idof] += flux[idof]; } } } else if (cells->is_local[l]) { - // Perform computation for a boundary edge PetscBool bnd_cell_order_flipped = PETSC_FALSE; @@ -949,12 +932,11 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { if (cells->centroid[l].X[0] > edges->centroid[iedge].X[0]) bnd_cell_order_flipped = PETSC_TRUE; } - PetscReal hl = x_ptr[l*dof + 0]; + PetscReal hl = x_ptr[l * dof + 0]; if (!(hl < user->tiny_h)) { - - PetscReal hul = x_ptr[l*dof + 1]; - PetscReal hvl = x_ptr[l*dof + 2]; + PetscReal hul = x_ptr[l * dof + 1]; + PetscReal hvl = x_ptr[l * dof + 2]; PetscReal ul, vl; PetscCall(GetVelocityFromMomentum(user->tiny_h, hl, hul, hvl, &ul, &vl)); @@ -962,28 +944,34 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { PetscReal hr, ur, vr; hr = hl; if (is_edge_vertical) { - ur = ul; + ur = ul; vr = -vl; } else { ur = -ul; - vr = vl; + vr = vl; } if (bnd_cell_order_flipped) { PetscReal tmp; - tmp = hl; hl = hr; hr = tmp; - tmp = ul; ul = ur; ur = tmp; - tmp = vl; vl = vr; vr = tmp; + tmp = hl; + hl = hr; + hr = tmp; + tmp = ul; + ul = ur; + ur = tmp; + tmp = vl; + vl = vr; + vr = tmp; } PetscReal flux[3], amax; PetscCall(solver(hl, hr, ul, ur, vl, vr, sn, cn, flux, &amax)); - for (PetscInt idof=0; idofNt * user->dt; - TS ts; + TS ts; PetscCall(TSCreate(user->comm, &ts)); PetscCall(TSSetProblemType(ts, TS_NONLINEAR)); PetscCall(TSSetType(ts, TSEULER)); @@ -1072,7 +1059,7 @@ int main(int argc, char **argv) { PetscCall(TSSetTimeStep(ts, user->dt)); PetscCall(TSSetFromOptions(ts)); - PetscCall(TSSolve(ts,X)); + PetscCall(TSSolve(ts, X)); if (user->save) { char fname[PETSC_MAX_PATH_LEN]; From 429a41f2efd9f64784dffabaa8136ec634732261 Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Wed, 9 Nov 2022 14:03:18 -0800 Subject: [PATCH 08/22] Adds the reflecting wall --- swe_roe/ex2.c | 183 ++++++++++++++++++++++++++++++++++++++++++++------ 1 file changed, 163 insertions(+), 20 deletions(-) diff --git a/swe_roe/ex2.c b/swe_roe/ex2.c index 58884d6..ed0b45f 100644 --- a/swe_roe/ex2.c +++ b/swe_roe/ex2.c @@ -120,12 +120,12 @@ struct _n_User { char filename[PETSC_MAX_PATH_LEN]; DM dm; PetscInt Nt, Nx, Ny; - PetscReal dt, hx, hy; + PetscReal dt, dx, dy; PetscReal Lx, Ly; PetscReal hu, hd; PetscReal tiny_h; PetscInt dof, rank, comm_size; - Vec B; + Vec B, localB; Vec localX; PetscBool debug, save, add_building; PetscInt tstep; @@ -167,8 +167,8 @@ static PetscErrorCode ProcessOptions(MPI_Comm comm, User user) { user->comm = comm; user->Nx = 4; user->Ny = 5; - user->hx = 1.0; - user->hy = 1.0; + 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; @@ -184,8 +184,8 @@ static PetscErrorCode ProcessOptions(MPI_Comm comm, User user) { 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("-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)); @@ -605,8 +605,6 @@ static PetscErrorCode SaveNaturalCellIDs(DM dm, RDyCell *cells, PetscInt rank) { // Add entries in the natural vector PetscScalar *entries; - printf("rank = %d; num_fields = %d; natural_size = %d %d; offset = %d\n", rank, num_fields, natural_size, cum_natural_size, - cum_natural_size / num_fields - natural_size / num_fields); PetscCall(VecGetArray(natural, &entries)); for (PetscInt i = 0; i < natural_size; ++i) { if (i % num_fields == 0) { @@ -636,7 +634,6 @@ static PetscErrorCode SaveNaturalCellIDs(DM dm, RDyCell *cells, PetscInt rank) { PetscInt local_size; PetscCall(VecGetLocalSize(local, &local_size)); PetscCall(VecGetArray(local, &entries)); - printf("rank = %d; local_size = %d; \n", rank, local_size); for (PetscInt i = 0; i < local_size / num_fields; ++i) { cells->natural_id[i] = entries[i * num_fields]; } @@ -685,7 +682,7 @@ static PetscErrorCode CreateMesh(User user) { PetscCall(PopulateCellsFromDM(user->dm, &mesh_ptr->cells, &mesh_ptr->num_cells_local)); PetscCall(PopulateEdgesFromDM(user->dm, &mesh_ptr->edges)); PetscCall(PopulateVerticesFromDM(user->dm, &mesh_ptr->vertices)); - PetscCall(SaveNaturalCellIDs(user->dm, &mesh_ptr->cells, user->rank)); + //PetscCall(SaveNaturalCellIDs(user->dm, &mesh_ptr->cells, user->rank)); PetscFunctionReturn(0); } @@ -720,6 +717,77 @@ static PetscErrorCode SetInitialCondition(User user, Vec X) { PetscFunctionReturn(0); } +PetscErrorCode AddBuildings(User user) { + PetscFunctionBeginUser; + + PetscInt bu = 30 / user->dx; + PetscInt bd = 105 / user->dx; + PetscInt bl = 95 / user->dy; + PetscInt br = 105 / user->dy; + + PetscCall(VecZeroEntries(user->B)); + + PetscScalar *b_ptr; + PetscCall(VecGetArray(user->B, &b_ptr)); + + /* + + x represents the reflecting wall, + o represents the dam that will be broken suddenly. + hu is the upsatream water dpeth, and hd is the downstream water depth. + + +/|\ xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx + | x xxxx x + | x xxxx x +95[m] x xxxx x + | x xxxx x + | x xxxx x + | x xxxx x +\|/ x xxxx x +/|\ x hu o hd x + | x o x + | x o x +105[m] x o x + | /|\ x xxxx x + | | x xxxx x + | 30[m] x xxxx x + | | x xxxx x + | | x xxxx x +\|/ \|/ xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx + | | | | + |<----- 95[m] ---->| | | + |<------- 105[m] ----->| | + |<---------------- 200[m] -------------->| + */ + + RDyMesh *mesh = user->mesh; + RDyCell *cells = &mesh->cells; + + PetscInt nbnd_1 = 0, nbnd_2 = 0; + for (PetscInt icell=0; icellnum_cells_local; icell++) { + PetscReal xc = cells->centroid[icell].X[1]; + PetscReal yc = cells->centroid[icell].X[0]; + if (yc < bu && xc >= bl && xc < br) { + b_ptr[icell*3 + 0] = 1.0; + nbnd_1++; + } else if (yc >= bd && xc >= bl && xc < br) { + b_ptr[icell*3 + 0] = 1.0; + nbnd_2++; + } + } + + PetscCall(VecRestoreArray(user->B, &b_ptr)); + + PetscCall(DMGlobalToLocalBegin(user->dm, user->B, INSERT_VALUES, user->localB)); + PetscCall(DMGlobalToLocalEnd(user->dm, user->B, INSERT_VALUES, user->localB)); + + PetscCall(PetscPrintf(user->comm, "Building size: bu=%d,bd=%d,bl=%d,br=%d\n", bu, bd, bl, br)); + PetscCall(PetscPrintf(user->comm, "Buildings added sucessfully!\n")); + + PetscFunctionReturn(0); +} + /// @brief Computes flux based on Roe solver /// @param [in] hl Height left of the edge /// @param [in] hr Height right of the edge @@ -754,7 +822,6 @@ PetscErrorCode solver(PetscReal hl, PetscReal hr, PetscReal ul, PetscReal ur, Pe PetscReal dv = vr - vl; PetscReal dupar = -du * sn + dv * cn; PetscReal duperp = du * cn + dv * sn; - // printf("cn = %f; sn = %f\n",cn,sn); PetscReal dW[3]; dW[0] = 0.5 * (dh - hhat * duperp / chat); @@ -870,9 +937,10 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { PetscCall(VecZeroEntries(F)); // Get pointers to vector data - PetscScalar *x_ptr, *f_ptr; + PetscScalar *x_ptr, *f_ptr, *b_ptr; PetscCall(VecGetArray(user->localX, &x_ptr)); PetscCall(VecGetArray(F, &f_ptr)); + PetscCall(VecGetArray(user->localB, &b_ptr)); PetscInt dof = 3; for (PetscInt iedge = 0; iedge < mesh->num_edges; iedge++) { @@ -900,28 +968,88 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { PetscReal hl = x_ptr[l * dof + 0]; PetscReal hr = x_ptr[r * dof + 0]; + PetscReal bl = b_ptr[l * dof + 0]; + PetscReal br = b_ptr[r * dof + 0]; - if (!(hr < user->tiny_h && hl < user->tiny_h)) { - PetscReal hul = x_ptr[l * dof + 1]; - PetscReal hvl = x_ptr[l * dof + 2]; + if (bl == 0 && br == 0) { + // Both, left and right cells are not boundary walls + if (!(hr < user->tiny_h && hl < user->tiny_h)) { + PetscReal hul = x_ptr[l * dof + 1]; + PetscReal hvl = x_ptr[l * dof + 2]; + PetscReal hur = x_ptr[r * dof + 1]; + PetscReal hvr = x_ptr[r * dof + 2]; + + PetscReal ur, vr, ul, vl; + + PetscCall(GetVelocityFromMomentum(user->tiny_h, hr, hur, hvr, &ur, &vr)); + PetscCall(GetVelocityFromMomentum(user->tiny_h, hl, hul, hvl, &ul, &vl)); + + PetscReal flux[3], amax; + PetscCall(solver(hl, hr, ul, ur, vl, vr, sn, cn, flux, &amax)); + + for (PetscInt idof = 0; idof < dof; idof++) { + if (cells->is_local[l]) f_ptr[l * dof + idof] -= flux[idof]; + if (cells->is_local[r]) f_ptr[r * dof + idof] += flux[idof]; + } + + } + + } else if (bl == 1 && br == 0) { + + // Left cell is a boundary wall and right cell is an internal cell + + PetscReal hr = x_ptr[r * dof + 0]; PetscReal hur = x_ptr[r * dof + 1]; PetscReal hvr = x_ptr[r * dof + 2]; - PetscReal ur, vr, ul, vl; - + PetscReal ur, vr; PetscCall(GetVelocityFromMomentum(user->tiny_h, hr, hur, hvr, &ur, &vr)); - PetscCall(GetVelocityFromMomentum(user->tiny_h, hl, hul, hvl, &ul, &vl)); + + PetscReal hl = hr; + PetscReal ul, vl; + if (is_edge_vertical) { + ul = ur; + vl = -vr; + } else { + ul = -ur; + vl = vr; + } PetscReal flux[3], amax; PetscCall(solver(hl, hr, ul, ur, vl, vr, sn, cn, flux, &amax)); + for (PetscInt idof = 0; idof < dof; idof++) { + if (cells->is_local[r]) f_ptr[r * dof + idof] += flux[idof]; + } + + } else if (bl == 0 && br == 1 ) { + // Left cell is an internal cell and right cell is a boundary wall + + PetscReal hl = x_ptr[l * dof + 0]; + PetscReal hul = x_ptr[l * dof + 1]; + PetscReal hvl = x_ptr[l * dof + 2]; + + PetscReal ul, vl; + PetscCall(GetVelocityFromMomentum(user->tiny_h, hl, hul, hvl, &ul, &vl)); + + PetscReal hr = hl; + PetscReal ur, vr; + if (is_edge_vertical) { + ur = ul; + vr = -vl; + } else { + ur = -ul; + vr = vl; + } + PetscReal flux[3], amax; + PetscCall(solver(hl, hr, ul, ur, vl, vr, sn, cn, flux, &amax)); for (PetscInt idof = 0; idof < dof; idof++) { if (cells->is_local[l]) f_ptr[l * dof + idof] -= flux[idof]; - if (cells->is_local[r]) f_ptr[r * dof + idof] += flux[idof]; } + } - } else if (cells->is_local[l]) { + } else if (cells->is_local[l] && b_ptr[l * dof + 0] == 0) { // Perform computation for a boundary edge PetscBool bnd_cell_order_flipped = PETSC_FALSE; @@ -981,6 +1109,7 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { // Restore vectors PetscCall(VecRestoreArray(user->localX, &x_ptr)); PetscCall(VecRestoreArray(F, &f_ptr)); + PetscCall(VecRestoreArray(user->localB, &b_ptr)); if (user->save) { char fname[PETSC_MAX_PATH_LEN]; @@ -994,6 +1123,12 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { PetscCall(PetscViewerBinaryOpen(user->comm, fname, FILE_MODE_WRITE, &viewer)); PetscCall(VecView(F, viewer)); PetscCall(PetscViewerDestroy(&viewer)); + + sprintf(fname, "outputs/ex2_b.dat"); + PetscCall(PetscViewerBinaryOpen(user->comm, fname, FILE_MODE_WRITE, &viewer)); + PetscCall(VecView(user->B, viewer)); + PetscCall(PetscViewerDestroy(&viewer)); + } PetscFunctionReturn(0); @@ -1021,6 +1156,7 @@ int main(int argc, char **argv) { PetscCall(VecDuplicate(X, &R)); VecViewFromOptions(X, NULL, "-vec_view"); PetscCall(DMCreateLocalVector(user->dm, &user->localX)); // size = dof * number of cells + PetscCall(VecDuplicate(user->localX, &user->localB)); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * * Create the mesh @@ -1043,6 +1179,13 @@ int main(int argc, char **argv) { PetscCall(PetscViewerDestroy(&viewer)); } + /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * + * Add buildings + * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ + if (user->add_building) { + PetscCall(AddBuildings(user)); + } + /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * * Create timestepping solver context * - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - */ From 999b203a491ab0703ee1ab0f7335d6f464f0e8d1 Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Wed, 9 Nov 2022 14:12:08 -0800 Subject: [PATCH 09/22] Accounts for edge length and cell area --- swe_roe/ex2.c | 20 ++++++++++++++------ 1 file changed, 14 insertions(+), 6 deletions(-) diff --git a/swe_roe/ex2.c b/swe_roe/ex2.c index ed0b45f..5239cf0 100644 --- a/swe_roe/ex2.c +++ b/swe_roe/ex2.c @@ -947,6 +947,7 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { PetscInt cellOffset = edges->cell_offset[iedge]; PetscInt l = edges->cell_ids[cellOffset]; PetscInt r = edges->cell_ids[cellOffset + 1]; + PetscReal edgeLen = edges->length[iedge]; PetscReal sn, cn; @@ -978,6 +979,8 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { PetscReal hvl = x_ptr[l * dof + 2]; PetscReal hur = x_ptr[r * dof + 1]; PetscReal hvr = x_ptr[r * dof + 2]; + PetscReal areal = cells->area[l]; + PetscReal arear = cells->area[r]; PetscReal ur, vr, ul, vl; @@ -988,8 +991,8 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { PetscCall(solver(hl, hr, ul, ur, vl, vr, sn, cn, flux, &amax)); for (PetscInt idof = 0; idof < dof; idof++) { - if (cells->is_local[l]) f_ptr[l * dof + idof] -= flux[idof]; - if (cells->is_local[r]) f_ptr[r * dof + idof] += flux[idof]; + if (cells->is_local[l]) f_ptr[l * dof + idof] -= flux[idof]*edgeLen/areal; + if (cells->is_local[r]) f_ptr[r * dof + idof] += flux[idof]*edgeLen/arear; } } @@ -1017,8 +1020,10 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { PetscReal flux[3], amax; PetscCall(solver(hl, hr, ul, ur, vl, vr, sn, cn, flux, &amax)); + + PetscReal arear = cells->area[r]; for (PetscInt idof = 0; idof < dof; idof++) { - if (cells->is_local[r]) f_ptr[r * dof + idof] += flux[idof]; + if (cells->is_local[r]) f_ptr[r * dof + idof] += flux[idof]*edgeLen/arear; } } else if (bl == 0 && br == 1 ) { @@ -1043,8 +1048,10 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { PetscReal flux[3], amax; PetscCall(solver(hl, hr, ul, ur, vl, vr, sn, cn, flux, &amax)); + + PetscReal areal = cells->area[l]; for (PetscInt idof = 0; idof < dof; idof++) { - if (cells->is_local[l]) f_ptr[l * dof + idof] -= flux[idof]; + if (cells->is_local[l]) f_ptr[l * dof + idof] -= flux[idof]*edgeLen/areal; } } @@ -1095,11 +1102,12 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { PetscReal flux[3], amax; PetscCall(solver(hl, hr, ul, ur, vl, vr, sn, cn, flux, &amax)); + PetscReal areal = cells->area[l]; for (PetscInt idof = 0; idof < dof; idof++) { if (!bnd_cell_order_flipped) { - f_ptr[l * dof + idof] -= flux[idof]; + f_ptr[l * dof + idof] -= flux[idof]*edgeLen/areal; } else { - f_ptr[l * dof + idof] += flux[idof]; + f_ptr[l * dof + idof] += flux[idof]*edgeLen/areal; } } } From 9dd9f78af04ecec9c153d6cdbb319519d1cc818b Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Wed, 9 Nov 2022 14:12:40 -0800 Subject: [PATCH 10/22] Fix formatting --- swe_roe/ex2.c | 56 ++++++++++++++++++++++++--------------------------- 1 file changed, 26 insertions(+), 30 deletions(-) diff --git a/swe_roe/ex2.c b/swe_roe/ex2.c index 5239cf0..7bee519 100644 --- a/swe_roe/ex2.c +++ b/swe_roe/ex2.c @@ -682,7 +682,7 @@ static PetscErrorCode CreateMesh(User user) { PetscCall(PopulateCellsFromDM(user->dm, &mesh_ptr->cells, &mesh_ptr->num_cells_local)); PetscCall(PopulateEdgesFromDM(user->dm, &mesh_ptr->edges)); PetscCall(PopulateVerticesFromDM(user->dm, &mesh_ptr->vertices)); - //PetscCall(SaveNaturalCellIDs(user->dm, &mesh_ptr->cells, user->rank)); + // PetscCall(SaveNaturalCellIDs(user->dm, &mesh_ptr->cells, user->rank)); PetscFunctionReturn(0); } @@ -761,18 +761,18 @@ PetscErrorCode AddBuildings(User user) { |<---------------- 200[m] -------------->| */ - RDyMesh *mesh = user->mesh; + RDyMesh *mesh = user->mesh; RDyCell *cells = &mesh->cells; PetscInt nbnd_1 = 0, nbnd_2 = 0; - for (PetscInt icell=0; icellnum_cells_local; icell++) { + for (PetscInt icell = 0; icell < mesh->num_cells_local; icell++) { PetscReal xc = cells->centroid[icell].X[1]; PetscReal yc = cells->centroid[icell].X[0]; if (yc < bu && xc >= bl && xc < br) { - b_ptr[icell*3 + 0] = 1.0; + b_ptr[icell * 3 + 0] = 1.0; nbnd_1++; } else if (yc >= bd && xc >= bl && xc < br) { - b_ptr[icell*3 + 0] = 1.0; + b_ptr[icell * 3 + 0] = 1.0; nbnd_2++; } } @@ -944,10 +944,10 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { PetscInt dof = 3; for (PetscInt iedge = 0; iedge < mesh->num_edges; iedge++) { - PetscInt cellOffset = edges->cell_offset[iedge]; - PetscInt l = edges->cell_ids[cellOffset]; - PetscInt r = edges->cell_ids[cellOffset + 1]; - PetscReal edgeLen = edges->length[iedge]; + PetscInt cellOffset = edges->cell_offset[iedge]; + PetscInt l = edges->cell_ids[cellOffset]; + PetscInt r = edges->cell_ids[cellOffset + 1]; + PetscReal edgeLen = edges->length[iedge]; PetscReal sn, cn; @@ -975,10 +975,10 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { if (bl == 0 && br == 0) { // Both, left and right cells are not boundary walls if (!(hr < user->tiny_h && hl < user->tiny_h)) { - PetscReal hul = x_ptr[l * dof + 1]; - PetscReal hvl = x_ptr[l * dof + 2]; - PetscReal hur = x_ptr[r * dof + 1]; - PetscReal hvr = x_ptr[r * dof + 2]; + PetscReal hul = x_ptr[l * dof + 1]; + PetscReal hvl = x_ptr[l * dof + 2]; + PetscReal hur = x_ptr[r * dof + 1]; + PetscReal hvr = x_ptr[r * dof + 2]; PetscReal areal = cells->area[l]; PetscReal arear = cells->area[r]; @@ -991,17 +991,15 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { PetscCall(solver(hl, hr, ul, ur, vl, vr, sn, cn, flux, &amax)); for (PetscInt idof = 0; idof < dof; idof++) { - if (cells->is_local[l]) f_ptr[l * dof + idof] -= flux[idof]*edgeLen/areal; - if (cells->is_local[r]) f_ptr[r * dof + idof] += flux[idof]*edgeLen/arear; + if (cells->is_local[l]) f_ptr[l * dof + idof] -= flux[idof] * edgeLen / areal; + if (cells->is_local[r]) f_ptr[r * dof + idof] += flux[idof] * edgeLen / arear; } - } } else if (bl == 1 && br == 0) { - // Left cell is a boundary wall and right cell is an internal cell - PetscReal hr = x_ptr[r * dof + 0]; + PetscReal hr = x_ptr[r * dof + 0]; PetscReal hur = x_ptr[r * dof + 1]; PetscReal hvr = x_ptr[r * dof + 2]; @@ -1011,11 +1009,11 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { PetscReal hl = hr; PetscReal ul, vl; if (is_edge_vertical) { - ul = ur; + ul = ur; vl = -vr; } else { ul = -ur; - vl = vr; + vl = vr; } PetscReal flux[3], amax; @@ -1023,13 +1021,13 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { PetscReal arear = cells->area[r]; for (PetscInt idof = 0; idof < dof; idof++) { - if (cells->is_local[r]) f_ptr[r * dof + idof] += flux[idof]*edgeLen/arear; + if (cells->is_local[r]) f_ptr[r * dof + idof] += flux[idof] * edgeLen / arear; } - } else if (bl == 0 && br == 1 ) { + } else if (bl == 0 && br == 1) { // Left cell is an internal cell and right cell is a boundary wall - PetscReal hl = x_ptr[l * dof + 0]; + PetscReal hl = x_ptr[l * dof + 0]; PetscReal hul = x_ptr[l * dof + 1]; PetscReal hvl = x_ptr[l * dof + 2]; @@ -1039,11 +1037,11 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { PetscReal hr = hl; PetscReal ur, vr; if (is_edge_vertical) { - ur = ul; + ur = ul; vr = -vl; } else { ur = -ul; - vr = vl; + vr = vl; } PetscReal flux[3], amax; @@ -1051,9 +1049,8 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { PetscReal areal = cells->area[l]; for (PetscInt idof = 0; idof < dof; idof++) { - if (cells->is_local[l]) f_ptr[l * dof + idof] -= flux[idof]*edgeLen/areal; + if (cells->is_local[l]) f_ptr[l * dof + idof] -= flux[idof] * edgeLen / areal; } - } } else if (cells->is_local[l] && b_ptr[l * dof + 0] == 0) { @@ -1105,9 +1102,9 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { PetscReal areal = cells->area[l]; for (PetscInt idof = 0; idof < dof; idof++) { if (!bnd_cell_order_flipped) { - f_ptr[l * dof + idof] -= flux[idof]*edgeLen/areal; + f_ptr[l * dof + idof] -= flux[idof] * edgeLen / areal; } else { - f_ptr[l * dof + idof] += flux[idof]*edgeLen/areal; + f_ptr[l * dof + idof] += flux[idof] * edgeLen / areal; } } } @@ -1136,7 +1133,6 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { PetscCall(PetscViewerBinaryOpen(user->comm, fname, FILE_MODE_WRITE, &viewer)); PetscCall(VecView(user->B, viewer)); PetscCall(PetscViewerDestroy(&viewer)); - } PetscFunctionReturn(0); From ea8d1f86bd8546ccc138834898dbd7e1d260d481 Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Wed, 9 Nov 2022 14:38:03 -0800 Subject: [PATCH 11/22] Changes the output filename --- swe_roe/ex1.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/swe_roe/ex1.c b/swe_roe/ex1.c index 2ae6bde..d949108 100644 --- a/swe_roe/ex1.c +++ b/swe_roe/ex1.c @@ -205,7 +205,7 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { 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); + sprintf(fname, "outputs/ex1_Nx_%d_Ny_%d_dt_%f_%d.dat", user->Nx, user->Ny, user->dt, user->tstep-1); PetscViewer viewer; PetscCall(PetscViewerBinaryOpen(user->comm, fname, FILE_MODE_WRITE, &viewer)); @@ -597,7 +597,7 @@ int main(int argc, char **argv) { 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); + PetscPrintf(user->comm, "Max simulation time is %f; \n", user->max_time); /* - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - * * Initialize DMDA @@ -698,7 +698,7 @@ int main(int argc, char **argv) { 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); + sprintf(fname, "outputs/ex1_Nx_%d_Ny_%d_dt_%f_%d.dat", user->Nx, user->Ny, user->dt, user->tstep); PetscViewer viewer; PetscCall(PetscViewerBinaryOpen(user->comm, fname, FILE_MODE_WRITE, &viewer)); From dfb325f3774f52d7a54e0a15037537edce9bc970 Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Wed, 9 Nov 2022 14:43:10 -0800 Subject: [PATCH 12/22] Updates the MATLAB script to visualize the output --- swe_roe/ex1_show.m | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/swe_roe/ex1_show.m b/swe_roe/ex1_show.m index 904070a..a80a1ed 100644 --- a/swe_roe/ex1_show.m +++ b/swe_roe/ex1_show.m @@ -1,11 +1,13 @@ clear;close all;clc -% User need to provide petsc directory -addpath('/Users/xudo627/developments/petsc/share/petsc/matlab/'); +if ~exist('PetscBinaryRead') + error(['Please add PETSc MATLAB files before calling this script via: ' char(10) ... + 'addpath /share/petsc/matlab/']) +end Lx = 200; % [m] Ly = 200; % [m] -Nt = length(dir('outputs/ex1*.dat')) - 1; +Nt = length(dir('outputs/ex1_Nx*.dat')) - 1; dof = 3; file_IC = dir('outputs/ex1_*_IC.dat'); @@ -43,7 +45,7 @@ imAlpha = ones(size(h0)); 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 +for i = 0 : Nt-1 file = dir(['outputs/ex1_*_' num2str(i) '.dat']); data = PetscBinaryRead(fullfile(file(1).folder,file(1).name)); data = reshape(data, [dof, length(data)/dof]); From 295ec128a1d8b461c2f4e755e3926620291126da Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Wed, 9 Nov 2022 14:53:52 -0800 Subject: [PATCH 13/22] Generalizes MATLAB plotting script --- swe_roe/{ex1_show.m => show_results.m} | 26 ++++++++++++++++++++------ 1 file changed, 20 insertions(+), 6 deletions(-) rename swe_roe/{ex1_show.m => show_results.m} (77%) diff --git a/swe_roe/ex1_show.m b/swe_roe/show_results.m similarity index 77% rename from swe_roe/ex1_show.m rename to swe_roe/show_results.m index a80a1ed..a52acd1 100644 --- a/swe_roe/ex1_show.m +++ b/swe_roe/show_results.m @@ -1,16 +1,28 @@ -clear;close all;clc +function show_results(example_number) +% Plots output of SWE from examples in this directory +% +% Examples +% show_results(1) + +close all if ~exist('PetscBinaryRead') error(['Please add PETSc MATLAB files before calling this script via: ' char(10) ... 'addpath /share/petsc/matlab/']) end +switch example_number + case {1,2} + otherwise + error('Invalid example_number. It can only be 1 or 2'); +end + Lx = 200; % [m] Ly = 200; % [m] -Nt = length(dir('outputs/ex1_Nx*.dat')) - 1; +Nt = length(dir(sprintf('outputs/ex%d_Nx*.dat',example_number))) - 1; dof = 3; -file_IC = dir('outputs/ex1_*_IC.dat'); +file_IC = dir(sprintf('outputs/ex%d_*_IC.dat',example_number)); if length(file_IC) > 1 error('Check your outputs, there may be two simulations!'); end @@ -46,7 +58,7 @@ imAlpha(1:30/dx,95/dy+1:105/dy) = 0; imAlpha(105/dy+1:200/dy,95/dy+1:105/dy) = 0; for i = 0 : Nt-1 - file = dir(['outputs/ex1_*_' num2str(i) '.dat']); + file = dir(['outputs/ex' num2str(example_number) '_*_' num2str(i) '.dat']); data = PetscBinaryRead(fullfile(file(1).folder,file(1).name)); data = reshape(data, [dof, length(data)/dof]); h = data(1,:); h = reshape(h,[Nx Ny]); @@ -54,7 +66,7 @@ v = data(3,:); v = reshape(v,[Nx Ny]); pause(0.01); imagesc(h,'AlphaData',imAlpha); colorbar; caxis([0 10]); - title(['t = ' num2str(i*dt) 's'],'FontSize',15); + title([ 'ex' num2str(example_number) ': t = ' num2str(i*dt) 's'],'FontSize',15); end h(1:30/dx+1,95/dy:105/dy+1) = NaN; @@ -71,6 +83,7 @@ colormap(jet); contour3(x,y,flipud(h),20, 'k-', 'LineWidth',1); view(30,45); +title(['ex' num2str(example_number)]) subplot(1,2,2); quiver(x,y,flipud(v),flipud(u)); hold on; @@ -78,4 +91,5 @@ fill([95 95 106 106 95],[170 200 200 170 170],[0.5 0.5 0.5],'EdgeColor','none'); xlim([1 200]);ylim([1 200]); contour(x,y,flipud(h),20,'LineWidth',1); -colorbar; colormap(jet); \ No newline at end of file +colorbar; colormap(jet); +title(['ex' num2str(example_number)]) From a967f90e7b3b57be72f932b51aeba0fe2d207c59 Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Wed, 9 Nov 2022 14:54:33 -0800 Subject: [PATCH 14/22] Minor modification to output from ex2 --- swe_roe/ex2.c | 11 +++-------- 1 file changed, 3 insertions(+), 8 deletions(-) diff --git a/swe_roe/ex2.c b/swe_roe/ex2.c index 7bee519..820a206 100644 --- a/swe_roe/ex2.c +++ b/swe_roe/ex2.c @@ -1118,21 +1118,16 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { if (user->save) { char fname[PETSC_MAX_PATH_LEN]; - sprintf(fname, "outputs/ex2_Nx_%d_Ny_%d_dt_%f_%d.dat", user->Nx, user->Ny, user->dt, user->tstep); + sprintf(fname, "outputs/ex2_Nx_%d_Ny_%d_dt_%f_%d.dat", user->Nx, user->Ny, user->dt, user->tstep-1); PetscViewer viewer; PetscCall(PetscViewerBinaryOpen(user->comm, fname, FILE_MODE_WRITE, &viewer)); PetscCall(VecView(X, viewer)); PetscCall(PetscViewerDestroy(&viewer)); - sprintf(fname, "outputs/ex2_flux_Nx_%d_Ny_%d_dt_%f_%d.dat", user->Nx, user->Ny, user->dt, user->tstep); + sprintf(fname, "outputs/ex2_flux_Nx_%d_Ny_%d_dt_%f_%d.dat", user->Nx, user->Ny, user->dt, user->tstep-1); PetscCall(PetscViewerBinaryOpen(user->comm, fname, FILE_MODE_WRITE, &viewer)); PetscCall(VecView(F, viewer)); PetscCall(PetscViewerDestroy(&viewer)); - - sprintf(fname, "outputs/ex2_b.dat"); - PetscCall(PetscViewerBinaryOpen(user->comm, fname, FILE_MODE_WRITE, &viewer)); - PetscCall(VecView(user->B, viewer)); - PetscCall(PetscViewerDestroy(&viewer)); } PetscFunctionReturn(0); @@ -1210,7 +1205,7 @@ int main(int argc, char **argv) { if (user->save) { char fname[PETSC_MAX_PATH_LEN]; - sprintf(fname, "outputs/ex2_Nx_%d_Ny_%d_dt_%f_final.dat", user->Nx, user->Ny, user->dt); + 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)); From 0d4916471bd3640725ca81442590d6a4961d6725 Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Wed, 9 Nov 2022 14:54:58 -0800 Subject: [PATCH 15/22] Fixes formatting --- swe_roe/ex1.c | 2 +- swe_roe/ex2.c | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/swe_roe/ex1.c b/swe_roe/ex1.c index d949108..d8fe416 100644 --- a/swe_roe/ex1.c +++ b/swe_roe/ex1.c @@ -205,7 +205,7 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { 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-1); + sprintf(fname, "outputs/ex1_Nx_%d_Ny_%d_dt_%f_%d.dat", user->Nx, user->Ny, user->dt, user->tstep - 1); PetscViewer viewer; PetscCall(PetscViewerBinaryOpen(user->comm, fname, FILE_MODE_WRITE, &viewer)); diff --git a/swe_roe/ex2.c b/swe_roe/ex2.c index 820a206..0f49760 100644 --- a/swe_roe/ex2.c +++ b/swe_roe/ex2.c @@ -1118,13 +1118,13 @@ PetscErrorCode RHSFunction(TS ts, PetscReal t, Vec X, Vec F, void *ptr) { if (user->save) { char fname[PETSC_MAX_PATH_LEN]; - sprintf(fname, "outputs/ex2_Nx_%d_Ny_%d_dt_%f_%d.dat", user->Nx, user->Ny, user->dt, user->tstep-1); + sprintf(fname, "outputs/ex2_Nx_%d_Ny_%d_dt_%f_%d.dat", user->Nx, user->Ny, user->dt, user->tstep - 1); PetscViewer viewer; PetscCall(PetscViewerBinaryOpen(user->comm, fname, FILE_MODE_WRITE, &viewer)); PetscCall(VecView(X, viewer)); PetscCall(PetscViewerDestroy(&viewer)); - sprintf(fname, "outputs/ex2_flux_Nx_%d_Ny_%d_dt_%f_%d.dat", user->Nx, user->Ny, user->dt, user->tstep-1); + sprintf(fname, "outputs/ex2_flux_Nx_%d_Ny_%d_dt_%f_%d.dat", user->Nx, user->Ny, user->dt, user->tstep - 1); PetscCall(PetscViewerBinaryOpen(user->comm, fname, FILE_MODE_WRITE, &viewer)); PetscCall(VecView(F, viewer)); PetscCall(PetscViewerDestroy(&viewer)); From bcd18ebcb3c92ccbc30f431227bc9e901afeb992 Mon Sep 17 00:00:00 2001 From: "Matthew G. Knepley" Date: Wed, 9 Nov 2022 19:34:56 -0500 Subject: [PATCH 16/22] I think the natural order works now --- swe_roe/ex2.c | 88 ++++++++++++++++++++++++++++++++++++++++----------- 1 file changed, 70 insertions(+), 18 deletions(-) diff --git a/swe_roe/ex2.c b/swe_roe/ex2.c index 0f49760..eb11c27 100644 --- a/swe_roe/ex2.c +++ b/swe_roe/ex2.c @@ -177,9 +177,7 @@ static PetscErrorCode ProcessOptions(MPI_Comm comm, User user) { MPI_Comm_size(user->comm, &user->comm_size); MPI_Comm_rank(user->comm, &user->rank); - PetscErrorCode ierr; - 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)); @@ -194,9 +192,8 @@ static PetscErrorCode ProcessOptions(MPI_Comm comm, User user) { PetscCall(PetscOptionsBool("-save", "save outputs", "", user->save, &user->save, NULL)); PetscCall(PetscOptionsString("-mesh_filename", "The mesh file", "ex2.c", user->filename, user->filename, PETSC_MAX_PATH_LEN, NULL)); } - ierr = PetscOptionsEnd(); + PetscOptionsEnd(); - PetscCall(ierr); assert(user->hu >= 0.); assert(user->hd >= 0.); @@ -597,38 +594,44 @@ static PetscErrorCode SaveNaturalCellIDs(DM dm, RDyCell *cells, PetscInt rank) { PetscCall(DMGetNumFields(dm, &num_fields)); // Create the natural vector - Vec natural; - PetscCall(DMCreateGlobalVector(dm, &natural)); - PetscInt natural_size = 0, cum_natural_size = 0; + Vec natural; + PetscInt natural_size = 0, natural_start; + PetscCall(DMPlexCreateNaturalVector(dm, &natural)); + PetscCall(PetscObjectSetName((PetscObject)natural, "Natural Vec")); PetscCall(VecGetLocalSize(natural, &natural_size)); - PetscCall(MPI_Scan(&natural_size, &cum_natural_size, 1, MPI_INT, MPI_SUM, PETSC_COMM_WORLD)); + PetscCall(VecGetOwnershipRange(natural, &natural_start, NULL)); // Add entries in the natural vector PetscScalar *entries; PetscCall(VecGetArray(natural, &entries)); for (PetscInt i = 0; i < natural_size; ++i) { if (i % num_fields == 0) { - entries[i] = i / num_fields + cum_natural_size / num_fields - natural_size / num_fields; + entries[i] = (natural_start + i) / num_fields; } else { - entries[i] = -1 - rank; + entries[i] = -(rank + 1); } } PetscCall(VecRestoreArray(natural, &entries)); - VecView(natural, PETSC_VIEWER_STDOUT_WORLD); + PetscCall(VecView(natural, PETSC_VIEWER_STDOUT_WORLD)); // Map natural IDs in global order Vec global; PetscCall(DMCreateGlobalVector(dm, &global)); + PetscCall(PetscObjectSetName((PetscObject)global, "Global Vec")); PetscCall(DMPlexNaturalToGlobalBegin(dm, natural, global)); PetscCall(DMPlexNaturalToGlobalEnd(dm, natural, global)); - VecView(global, PETSC_VIEWER_STDOUT_WORLD); + PetscCall(VecView(global, PETSC_VIEWER_STDOUT_WORLD)); // Map natural IDs in local order - Vec local; + Vec local; + PetscViewer selfviewer; PetscCall(DMCreateLocalVector(dm, &local)); + PetscCall(PetscObjectSetName((PetscObject)local, "Local Vec")); PetscCall(DMGlobalToLocalBegin(dm, global, INSERT_VALUES, local)); PetscCall(DMGlobalToLocalEnd(dm, global, INSERT_VALUES, local)); - VecView(local, PETSC_VIEWER_STDOUT_WORLD); + PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &selfviewer)); + PetscCall(VecView(local, selfviewer)); + PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &selfviewer)); // Save natural IDs PetscInt local_size; @@ -682,11 +685,54 @@ static PetscErrorCode CreateMesh(User user) { PetscCall(PopulateCellsFromDM(user->dm, &mesh_ptr->cells, &mesh_ptr->num_cells_local)); PetscCall(PopulateEdgesFromDM(user->dm, &mesh_ptr->edges)); PetscCall(PopulateVerticesFromDM(user->dm, &mesh_ptr->vertices)); - // PetscCall(SaveNaturalCellIDs(user->dm, &mesh_ptr->cells, user->rank)); + PetscCall(SaveNaturalCellIDs(user->dm, &mesh_ptr->cells, user->rank)); PetscFunctionReturn(0); } +static PetscErrorCode DestroyMesh(User user) { + PetscFunctionBegin; + PetscCall(PetscFree(user->mesh->cells.id)); + PetscCall(PetscFree(user->mesh->cells.global_id)); + PetscCall(PetscFree(user->mesh->cells.natural_id)); + PetscCall(PetscFree(user->mesh->cells.is_local)); + PetscCall(PetscFree(user->mesh->cells.num_vertices)); + PetscCall(PetscFree(user->mesh->cells.num_edges)); + PetscCall(PetscFree(user->mesh->cells.num_neighbors)); + PetscCall(PetscFree(user->mesh->cells.vertex_offset)); + PetscCall(PetscFree(user->mesh->cells.edge_offset)); + PetscCall(PetscFree(user->mesh->cells.neighbor_offset)); + PetscCall(PetscFree(user->mesh->cells.vertex_ids)); + PetscCall(PetscFree(user->mesh->cells.edge_ids)); + PetscCall(PetscFree(user->mesh->cells.neighbor_ids)); + PetscCall(PetscFree(user->mesh->cells.centroid)); + PetscCall(PetscFree(user->mesh->cells.area)); + PetscCall(PetscFree(user->mesh->edges.id)); + PetscCall(PetscFree(user->mesh->edges.global_id)); + PetscCall(PetscFree(user->mesh->edges.is_local)); + PetscCall(PetscFree(user->mesh->edges.num_cells)); + PetscCall(PetscFree(user->mesh->edges.vertex_ids)); + PetscCall(PetscFree(user->mesh->edges.cell_offset)); + PetscCall(PetscFree(user->mesh->edges.cell_ids)); + PetscCall(PetscFree(user->mesh->edges.is_internal)); + PetscCall(PetscFree(user->mesh->edges.normal)); + PetscCall(PetscFree(user->mesh->edges.centroid)); + PetscCall(PetscFree(user->mesh->edges.length)); + PetscCall(PetscFree(user->mesh->vertices.id)); + PetscCall(PetscFree(user->mesh->vertices.global_id)); + PetscCall(PetscFree(user->mesh->vertices.is_local)); + PetscCall(PetscFree(user->mesh->vertices.num_cells)); + PetscCall(PetscFree(user->mesh->vertices.num_edges)); + PetscCall(PetscFree(user->mesh->vertices.cell_offset)); + PetscCall(PetscFree(user->mesh->vertices.edge_offset)); + PetscCall(PetscFree(user->mesh->vertices.cell_ids)); + PetscCall(PetscFree(user->mesh->vertices.edge_ids)); + PetscCall(PetscFree(user->mesh->vertices.coordinate)); + PetscCall(PetscFree(user->mesh->nG2A)); + PetscCall(PetscFree(user->mesh)); + PetscFunctionReturn(0); +} + /// @brief Sets initial condition for [h, hu, hv] /// @param [in] user A User data structure /// @param [inout] X Vec for initial condition @@ -1153,7 +1199,7 @@ int main(int argc, char **argv) { PetscCall(DMCreateGlobalVector(user->dm, &X)); // size = dof * number of cells PetscCall(VecDuplicate(X, &user->B)); PetscCall(VecDuplicate(X, &R)); - VecViewFromOptions(X, NULL, "-vec_view"); + PetscCall(VecViewFromOptions(X, NULL, "-vec_view")); PetscCall(DMCreateLocalVector(user->dm, &user->localX)); // size = dof * number of cells PetscCall(VecDuplicate(user->localX, &user->localB)); @@ -1213,10 +1259,16 @@ int main(int argc, char **argv) { PetscCall(PetscViewerDestroy(&viewer)); } + PetscCall(TSDestroy(&ts)); PetscCall(VecDestroy(&X)); PetscCall(VecDestroy(&R)); + PetscCall(VecDestroy(&user->B)); + PetscCall(VecDestroy(&user->localB)); + PetscCall(VecDestroy(&user->localX)); + PetscCall(VecDestroy(&R)); + PetscCall(DestroyMesh(user)); PetscCall(DMDestroy(&user->dm)); - PetscCall(TSDestroy(&ts)); + PetscCall(PetscFree(user)); PetscCall(PetscFinalize()); From 17d07345e3a2547a4e8007c309685cb8271f5690 Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Thu, 10 Nov 2022 10:19:05 -0800 Subject: [PATCH 17/22] Minor fixes for PETScv3.18.1 --- swe_roe/ex2.c | 7 ++----- 1 file changed, 2 insertions(+), 5 deletions(-) diff --git a/swe_roe/ex2.c b/swe_roe/ex2.c index 0f49760..7149568 100644 --- a/swe_roe/ex2.c +++ b/swe_roe/ex2.c @@ -177,9 +177,7 @@ static PetscErrorCode ProcessOptions(MPI_Comm comm, User user) { MPI_Comm_size(user->comm, &user->comm_size); MPI_Comm_rank(user->comm, &user->rank); - PetscErrorCode ierr; - 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)); @@ -194,9 +192,8 @@ static PetscErrorCode ProcessOptions(MPI_Comm comm, User user) { PetscCall(PetscOptionsBool("-save", "save outputs", "", user->save, &user->save, NULL)); PetscCall(PetscOptionsString("-mesh_filename", "The mesh file", "ex2.c", user->filename, user->filename, PETSC_MAX_PATH_LEN, NULL)); } - ierr = PetscOptionsEnd(); + PetscOptionsEnd(); - PetscCall(ierr); assert(user->hu >= 0.); assert(user->hd >= 0.); From 926b8c0607e8857a2d9f2abb7f05954aa7c9be24 Mon Sep 17 00:00:00 2001 From: "Matthew G. Knepley" Date: Fri, 11 Nov 2022 08:57:44 -0500 Subject: [PATCH 18/22] Do not distribute by default --- swe_roe/ex2.c | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/swe_roe/ex2.c b/swe_roe/ex2.c index eb11c27..bce804b 100644 --- a/swe_roe/ex2.c +++ b/swe_roe/ex2.c @@ -227,9 +227,11 @@ static PetscErrorCode CreateDM(User user) { } else { DMPlexCreateFromFile(user->comm, user->filename, "ex2.c", PETSC_FALSE, &user->dm); } + PetscCall(DMPlexDistributeSetDefault(user->dm, PETSC_FALSE)); PetscObjectSetName((PetscObject)user->dm, "Mesh"); PetscCall(DMSetFromOptions(user->dm)); + PetscCall(DMViewFromOptions(user->dm, NULL, "-orig_dm_view")); // Determine the number of cells, edges, and vertices of the mesh PetscInt cStart, cEnd; @@ -278,6 +280,10 @@ static PetscErrorCode CreateDM(User user) { } PetscCall(DMViewFromOptions(user->dm, NULL, "-dm_view")); + PetscSF sf; + PetscCall(DMPlexGetGlobalToNaturalSF(user->dm, &sf)); + PetscCall(PetscObjectViewFromOptions((PetscObject)sf, NULL, "-nat_sf_view")); + PetscFunctionReturn(0); } From 9d2c08804a54654d88761a045b66e8ccf72f4e4a Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Tue, 15 Nov 2022 08:33:35 -0800 Subject: [PATCH 19/22] PetscSF code block --- swe_roe/ex2.c | 2 ++ 1 file changed, 2 insertions(+) diff --git a/swe_roe/ex2.c b/swe_roe/ex2.c index bce804b..a57f657 100644 --- a/swe_roe/ex2.c +++ b/swe_roe/ex2.c @@ -280,9 +280,11 @@ static PetscErrorCode CreateDM(User user) { } PetscCall(DMViewFromOptions(user->dm, NULL, "-dm_view")); +/* PetscSF sf; PetscCall(DMPlexGetGlobalToNaturalSF(user->dm, &sf)); PetscCall(PetscObjectViewFromOptions((PetscObject)sf, NULL, "-nat_sf_view")); +*/ PetscFunctionReturn(0); } From 107eeee2cfd629cb20424527f959d29c7bf37d03 Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Tue, 15 Nov 2022 08:34:50 -0800 Subject: [PATCH 20/22] Delete VecViews --- swe_roe/ex2.c | 3 --- 1 file changed, 3 deletions(-) diff --git a/swe_roe/ex2.c b/swe_roe/ex2.c index a57f657..3095593 100644 --- a/swe_roe/ex2.c +++ b/swe_roe/ex2.c @@ -620,7 +620,6 @@ static PetscErrorCode SaveNaturalCellIDs(DM dm, RDyCell *cells, PetscInt rank) { } } PetscCall(VecRestoreArray(natural, &entries)); - PetscCall(VecView(natural, PETSC_VIEWER_STDOUT_WORLD)); // Map natural IDs in global order Vec global; @@ -628,7 +627,6 @@ static PetscErrorCode SaveNaturalCellIDs(DM dm, RDyCell *cells, PetscInt rank) { PetscCall(PetscObjectSetName((PetscObject)global, "Global Vec")); PetscCall(DMPlexNaturalToGlobalBegin(dm, natural, global)); PetscCall(DMPlexNaturalToGlobalEnd(dm, natural, global)); - PetscCall(VecView(global, PETSC_VIEWER_STDOUT_WORLD)); // Map natural IDs in local order Vec local; @@ -638,7 +636,6 @@ static PetscErrorCode SaveNaturalCellIDs(DM dm, RDyCell *cells, PetscInt rank) { PetscCall(DMGlobalToLocalBegin(dm, global, INSERT_VALUES, local)); PetscCall(DMGlobalToLocalEnd(dm, global, INSERT_VALUES, local)); PetscCall(PetscViewerGetSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &selfviewer)); - PetscCall(VecView(local, selfviewer)); PetscCall(PetscViewerRestoreSubViewer(PETSC_VIEWER_STDOUT_WORLD, PETSC_COMM_SELF, &selfviewer)); // Save natural IDs From a6887fa0b6ba6686d6fc31389e574da2a1f5343f Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Tue, 15 Nov 2022 08:45:24 -0800 Subject: [PATCH 21/22] Address comments during review --- swe_roe/ex2.c | 8 +++----- swe_roe/show_results.m | 12 ++++++++++-- 2 files changed, 13 insertions(+), 7 deletions(-) diff --git a/swe_roe/ex2.c b/swe_roe/ex2.c index 3095593..1332182 100644 --- a/swe_roe/ex2.c +++ b/swe_roe/ex2.c @@ -197,13 +197,11 @@ static PetscErrorCode ProcessOptions(MPI_Comm comm, User user) { assert(user->hu >= 0.); assert(user->hd >= 0.); - user->Lx = user->Nx * 1.0; - user->Ly = user->Ny * 1.0; + user->Lx = user->Nx * user->dx; + user->Ly = user->Ny * user->dy; PetscReal max_time = user->Nt * user->dt; - if (!user->rank) { - PetscPrintf(user->comm, "Max simulation time is %f\n", max_time); - } + PetscPrintf(user->comm, "Max simulation time is %f\n", max_time); PetscFunctionReturn(0); } diff --git a/swe_roe/show_results.m b/swe_roe/show_results.m index a52acd1..0c5aa00 100644 --- a/swe_roe/show_results.m +++ b/swe_roe/show_results.m @@ -6,9 +6,17 @@ function show_results(example_number) close all + if ~exist('PetscBinaryRead') - error(['Please add PETSc MATLAB files before calling this script via: ' char(10) ... - 'addpath /share/petsc/matlab/']) + + % Let's try to get information about PETSc from environment variable + [~,PETSC_DIR] = system('echo $PETSC_DIR'); + if (exists(PETSC_DIR)) + addpath([PETSC_DIR(1:end-1) '/share/petsc/matlab/']) + else + error(['Please add PETSc MATLAB files before calling this script via: ' char(10) ... + 'addpath /share/petsc/matlab/']) + end end switch example_number From aa040fb9fbb20406c713bf1de05865069004c00a Mon Sep 17 00:00:00 2001 From: Gautam Bisht Date: Tue, 15 Nov 2022 08:52:43 -0800 Subject: [PATCH 22/22] Fixes formatting --- swe_roe/ex2.c | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/swe_roe/ex2.c b/swe_roe/ex2.c index 1332182..c6a3f6f 100644 --- a/swe_roe/ex2.c +++ b/swe_roe/ex2.c @@ -278,11 +278,11 @@ static PetscErrorCode CreateDM(User user) { } PetscCall(DMViewFromOptions(user->dm, NULL, "-dm_view")); -/* - PetscSF sf; - PetscCall(DMPlexGetGlobalToNaturalSF(user->dm, &sf)); - PetscCall(PetscObjectViewFromOptions((PetscObject)sf, NULL, "-nat_sf_view")); -*/ + /* + PetscSF sf; + PetscCall(DMPlexGetGlobalToNaturalSF(user->dm, &sf)); + PetscCall(PetscObjectViewFromOptions((PetscObject)sf, NULL, "-nat_sf_view")); + */ PetscFunctionReturn(0); }