diff --git a/docs/Learning/BoundaryConditions.md b/docs/Learning/BoundaryConditions.md new file mode 100644 index 000000000..cbf9c7926 --- /dev/null +++ b/docs/Learning/BoundaryConditions.md @@ -0,0 +1,260 @@ +# Boundary Condition System + +This page describes the internal design of SELF's extensible boundary condition (BC) system. It is intended for developers who want to understand the architecture, add new BC types to built-in models, or extend the system itself. + +For a practical guide to using boundary conditions in your own applications, see [Boundary Conditions (User Guide)](../Models/boundary-conditions.md). + +## Design Goals + +The BC system is designed to: + +- Allow each model to register only the BC types it supports +- Dispatch BC application through procedure pointers, avoiding monolithic `select case` blocks +- Pre-compute which mesh faces belong to each BC type, so BC routines loop only over relevant faces +- Support GPU-accelerated BC kernels with per-BC device arrays + +## Architecture Overview + +The system is built on three components: + +1. **`BoundaryCondition`** -- a node in a linked list that holds a BC identifier, a procedure pointer to the BC implementation, and arrays of element/side indices +2. **`BoundaryConditionList`** -- a doubly-linked list that manages registration, lookup, and iteration over `BoundaryCondition` nodes +3. **`SELF_bcMethod`** -- an abstract interface that all BC implementations must satisfy + +These are defined in `src/SELF_BoundaryConditions.f90`. + +### The `BoundaryCondition` Type + +```fortran +type BoundaryCondition + procedure(SELF_bcMethod), pointer :: bcMethod => null() + integer :: bcid + character(SELF_BCNAME_LENGTH) :: bcname + integer :: nBoundaries + integer, allocatable :: elements(:) + integer, allocatable :: sides(:) + type(c_ptr) :: elements_gpu = c_null_ptr + type(c_ptr) :: sides_gpu = c_null_ptr + type(BoundaryCondition), pointer :: next => null() + type(BoundaryCondition), pointer :: prev => null() +endtype +``` + +| Field | Purpose | +|-------|---------| +| `bcMethod` | Procedure pointer to the BC implementation | +| `bcid` | Integer constant identifying the BC type (e.g., `SELF_BC_PRESCRIBED`) | +| `bcname` | Human-readable name for diagnostics | +| `nBoundaries` | Number of boundary faces that carry this BC | +| `elements(:)` | Element indices for each boundary face | +| `sides(:)` | Local side indices for each boundary face | +| `elements_gpu`, `sides_gpu` | Device pointers used by GPU kernels | + +### The `SELF_bcMethod` Interface + +Every BC implementation must match this signature: + +```fortran +subroutine SELF_bcMethod(this, mymodel) + class(BoundaryCondition), intent(in) :: this + class(Model), intent(inout) :: mymodel +endsubroutine +``` + +The BC receives itself (providing access to `elements`, `sides`, and `nBoundaries`) and the model (providing access to solution data). Implementations use `select type` to downcast `mymodel` to the concrete model type. + +### The `BoundaryConditionList` Type + +Each `DGModel` carries two lists: + +- `hyperbolicBCs` -- for boundary conditions on the solution (used by `SetBoundaryCondition`) +- `parabolicBCs` -- for boundary conditions on the solution gradient (used by `SetGradientBoundaryCondition`) + +Key methods: + +| Method | Purpose | +|--------|---------| +| `Init()` | Initialize an empty list | +| `Free()` | Deallocate all nodes | +| `RegisterBoundaryCondition(bcid, bcname, bcfunc)` | Add a new BC or update an existing one | +| `GetBCForID(bcid)` | Return the node for a given `bcid`, or `null()` | +| `PopulateBoundaries(bcid, nBoundaries, elements, sides)` | Fill element/side arrays after mesh scanning | + +If `RegisterBoundaryCondition` is called with a `bcid` that is already registered, it updates the procedure pointer without creating a new node. This is how GPU model variants override CPU implementations. + +## BC Identifier Constants + +BC type identifiers are integer parameters defined in `src/SELF_Mesh.f90`: + +```fortran +! Conditions on the solution +integer, parameter :: SELF_BC_PRESCRIBED = 100 +integer, parameter :: SELF_BC_RADIATION = 101 +integer, parameter :: SELF_BC_NONORMALFLOW = 102 + +! Conditions on the solution gradients +integer, parameter :: SELF_BC_PRESCRIBED_STRESS = 200 +integer, parameter :: SELF_BC_NOSTRESS = 201 +``` + +These constants are used both when tagging mesh faces (in `sideInfo`) and when registering BCs. + +## Initialization Flow + +The BC system is initialized as part of model creation. The sequence is: + +``` +Model%Init(mesh, geometry) + | + +-- hyperbolicBCs%Init() + +-- parabolicBCs%Init() + +-- AdditionalInit() <-- subclass registers BCs here + +-- MapBoundaryConditions() <-- scans mesh, populates element/side arrays +``` + +### Step 1: Register BCs in `AdditionalInit` + +Each model subclass overrides `AdditionalInit` to register its supported BC types: + +```fortran +subroutine AdditionalInit_ECAdvection2D_t(this) + class(ECAdvection2D_t), intent(inout) :: this + procedure(SELF_bcMethod), pointer :: bcfunc + + bcfunc => hbc2d_NoNormalFlow_ECAdvection2D + call this%hyperbolicBCs%RegisterBoundaryCondition( & + SELF_BC_NONORMALFLOW, "no_normal_flow", bcfunc) +endsubroutine +``` + +### Step 2: Map Mesh Faces in `MapBoundaryConditions` + +After registration, `MapBoundaryConditions` scans the mesh `sideInfo` array. For each registered BC, it performs two passes: + +1. **Count** how many boundary faces carry that `bcid` +2. **Collect** the element and side indices into arrays + +These arrays are stored in the `BoundaryCondition` node via `PopulateBoundaries`. A boundary face is identified by `sideInfo(3,j,iEl) == 0` (no neighbor element) and `sideInfo(5,j,iEl) == bcid`. + +## Runtime Dispatch + +During time integration, `SetBoundaryCondition` iterates through the linked list and calls each registered BC: + +```fortran +subroutine SetBoundaryCondition(this) + class(DGModel2D_t), intent(inout) :: this + type(BoundaryCondition), pointer :: bc + procedure(SELF_bcMethod), pointer :: apply_bc + + bc => this%hyperbolicBCs%head + do while (associated(bc)) + apply_bc => bc%bcMethod + call apply_bc(bc, this) + bc => bc%next + enddo +endsubroutine +``` + +The same pattern applies to `SetGradientBoundaryCondition` for parabolic BCs. + +## GPU Acceleration + +GPU-enabled models follow a layered pattern: + +1. The CPU model class (e.g., `ECAdvection2D_t`) registers a Fortran BC implementation in its `AdditionalInit` +2. The GPU model class (e.g., `ECAdvection2D` in `src/gpu/`) extends the CPU class and: + - Calls the parent `AdditionalInit` to register the CPU version + - Re-registers the same `bcid` with a GPU wrapper function, which replaces the procedure pointer +3. During `Init`, the GPU class uploads `elements` and `sides` arrays to device memory (`elements_gpu`, `sides_gpu`) +4. During `Free`, the GPU class deallocates device arrays + +### GPU Wrapper Pattern + +A GPU wrapper is a Fortran subroutine matching `SELF_bcMethod` that calls a C/C++ kernel: + +```fortran +subroutine hbc2d_Mirror_ECAdvection2D_GPU_wrapper(bc, mymodel) + class(BoundaryCondition), intent(in) :: bc + class(Model), intent(inout) :: mymodel + + select type (m => mymodel) + class is (ECAdvection2D) + if (bc%nBoundaries > 0) then + call hbc2d_mirror_ecadvection2d_gpu( & + m%solution%extBoundary_gpu, & + m%solution%boundary_gpu, & + bc%elements_gpu, bc%sides_gpu, & + bc%nBoundaries, m%solution%interp%N, & + m%solution%nElem, m%solution%nvar) + endif + endselect +endsubroutine +``` + +The C++ kernel receives device pointers and iterates over the pre-filtered boundary face list: + +```cpp +__global__ void hbc2d_mirror_ecadvection2d_kernel( + real *extBoundary, real *boundary, + int *elements, int *sides, + int nBoundaries, int N, int nel, int nvar) +{ + // Thread indexing over DOFs, boundary faces, and variables + // elements[n] and sides[n] identify which face to process +} +``` + +### GPU Memory Lifecycle + +``` +Init_ECAdvection2D(mesh, geometry) + | + +-- Init_ECDGModel2D_t() (parent: registers CPU BCs, maps mesh) + +-- for each BC in hyperbolicBCs: + hipMalloc(elements_gpu) + hipMemcpy(elements -> elements_gpu) + hipMalloc(sides_gpu) + hipMemcpy(sides -> sides_gpu) + +Free_ECAdvection2D() + | + +-- for each BC in hyperbolicBCs: + hipFree(elements_gpu) + hipFree(sides_gpu) + +-- Free_ECDGModel2D_t() (parent: frees BC list nodes) +``` + +## Adding a New BC Type to a Built-in Model + +To add a new BC type (e.g., an inflow condition) to an existing model: + +1. **Define a BC ID** in `src/SELF_Mesh.f90` if one does not already exist: + + ```fortran + integer, parameter :: SELF_BC_INFLOW = 103 + ``` + +2. **Write the BC implementation** in the model's `_t` source file, matching the `SELF_bcMethod` interface + +3. **Register it** in the model's `AdditionalInit`: + + ```fortran + bcfunc => hbc2d_Inflow_MyModel + call this%hyperbolicBCs%RegisterBoundaryCondition( & + SELF_BC_INFLOW, "inflow", bcfunc) + ``` + +4. **For GPU models**, write a C++ kernel and Fortran wrapper, then re-register in the GPU class `AdditionalInit` + +5. **Tag mesh faces** with the new BC ID in the mesh setup (e.g., via `ResetBoundaryConditionType` or by setting `sideInfo(5,:,:)` appropriately) + +## Key Source Files + +| File | Contents | +|------|----------| +| `src/SELF_BoundaryConditions.f90` | `BoundaryCondition`, `BoundaryConditionList`, `SELF_bcMethod` interface | +| `src/SELF_Mesh.f90` | BC ID constants (`SELF_BC_PRESCRIBED`, etc.) | +| `src/SELF_DGModel{1D,2D,3D}_t.f90` | `MapBoundaryConditions`, `SetBoundaryCondition`, `SetGradientBoundaryCondition` | +| `src/SELF_Model.f90` | Base `AdditionalInit` / `AdditionalFree` stubs | +| `src/gpu/SELF_ECAdvection2D.f90` | Example GPU BC wrapper pattern | +| `src/gpu/SELF_ECAdvection2D.cpp` | Example GPU BC kernel | diff --git a/docs/Models/boundary-conditions.md b/docs/Models/boundary-conditions.md new file mode 100644 index 000000000..e17bf2a11 --- /dev/null +++ b/docs/Models/boundary-conditions.md @@ -0,0 +1,258 @@ +# Boundary Conditions + +SELF uses an extensible boundary condition system that lets you register custom boundary conditions for your models. This page explains how to use boundary conditions in practice. + +For details on the internal architecture, see [Boundary Condition System (Developer Guide)](../Learning/BoundaryConditions.md). + +## Overview + +Every DG model in SELF maintains two boundary condition lists: + +- **Hyperbolic BCs** -- conditions on the solution at domain boundaries +- **Parabolic BCs** -- conditions on the solution gradient at domain boundaries + +You populate these lists by overriding the `AdditionalInit` method in your model subclass. SELF then automatically determines which mesh faces match each BC and calls your BC routines during time integration. + +## Available BC Types + +SELF provides the following built-in BC identifiers: + +| Constant | Value | Category | Description | +|----------|-------|----------|-------------| +| `SELF_BC_PRESCRIBED` | 100 | Hyperbolic | Dirichlet-type: set the exterior state to a prescribed value | +| `SELF_BC_RADIATION` | 101 | Hyperbolic | Radiation / non-reflecting outflow | +| `SELF_BC_NONORMALFLOW` | 102 | Hyperbolic | Mirror / no-normal-flow (wall) | +| `SELF_BC_PRESCRIBED_STRESS` | 200 | Parabolic | Prescribed stress on the gradient | +| `SELF_BC_NOSTRESS` | 201 | Parabolic | Zero-stress (free-slip) | + +These constants are defined in `SELF_Mesh` and are used both for tagging mesh faces and for registering BC implementations. + +## Workflow + +Setting up boundary conditions involves three steps: + +1. **Tag mesh faces** with BC identifiers +2. **Write BC subroutines** that set exterior state values +3. **Register BCs** in your model's `AdditionalInit` + +### Step 1: Tag Mesh Faces + +After creating your mesh, assign BC identifiers to boundary faces. The simplest way is `ResetBoundaryConditionType`: + +```fortran +! Set left boundary to prescribed, right boundary to prescribed +call mesh%ResetBoundaryConditionType(SELF_BC_PRESCRIBED, SELF_BC_PRESCRIBED) +``` + +For 2D and 3D meshes, boundary face tagging is typically set via the mesh file or by modifying `sideInfo(5,:,:)` after mesh creation. + +### Step 2: Write BC Subroutines + +Each BC subroutine must match the `SELF_bcMethod` interface: + +```fortran +subroutine my_bc(bc, mymodel) + use SELF_BoundaryConditions + use SELF_Model + class(BoundaryCondition), intent(in) :: bc + class(Model), intent(inout) :: mymodel +endsubroutine +``` + +Inside the subroutine, use `select type` to access your model's data, then loop over `bc%nBoundaries` to set the exterior state on each boundary face: + +```fortran +subroutine hbc1d_Prescribed_mymodel(bc, mymodel) + class(BoundaryCondition), intent(in) :: bc + class(Model), intent(inout) :: mymodel + integer :: n, iEl, s + + select type (m => mymodel) + class is (my_model_type) + do n = 1, bc%nBoundaries + iEl = bc%elements(n) ! element index + s = bc%sides(n) ! local side index + ! Set the exterior boundary state + m%solution%extBoundary(s, iEl, 1:m%nvar) = ... + enddo + endselect +endsubroutine +``` + +!!! note + The `extBoundary` array holds the exterior (ghost) state used by the Riemann solver at domain boundaries. Setting `extBoundary = boundary` (the interior state) produces a mirror/no-normal-flow condition. Setting it to a specific value produces a Dirichlet condition. + +### Step 3: Register BCs in `AdditionalInit` + +Override `AdditionalInit` in your model subclass and register each BC: + +```fortran +subroutine AdditionalInit_mymodel(this) + class(my_model_type), intent(inout) :: this + procedure(SELF_bcMethod), pointer :: bcfunc + + ! Register a hyperbolic BC + bcfunc => hbc1d_Prescribed_mymodel + call this%hyperbolicBCs%RegisterBoundaryCondition( & + SELF_BC_PRESCRIBED, "prescribed", bcfunc) + + ! Register a parabolic BC (if needed) + bcfunc => pbc1d_Prescribed_mymodel + call this%parabolicBCs%RegisterBoundaryCondition( & + SELF_BC_PRESCRIBED, "prescribed", bcfunc) +endsubroutine +``` + +SELF calls `AdditionalInit` during model initialization, before scanning the mesh. After your BCs are registered, SELF automatically determines which mesh faces belong to each BC type and stores the element/side arrays. + +## Complete Example: Traveling Shock (Burgers 1D) + +This example shows how to extend the built-in `burgers1D` model with prescribed boundary conditions for a traveling shock solution. + +### Model Module + +```fortran +module burgers1d_shock_model + use self_Burgers1D + use SELF_BoundaryConditions + implicit none + + type, extends(burgers1D) :: burgers1d_shock + real(prec) :: ul = 1.0_prec ! Left state + real(prec) :: ur = 0.0_prec ! Right state + real(prec) :: x0 = 0.1_prec ! Initial shock position + contains + procedure :: AdditionalInit => AdditionalInit_burgers1d_shock + endtype + +contains + + subroutine AdditionalInit_burgers1d_shock(this) + class(burgers1d_shock), intent(inout) :: this + procedure(SELF_bcMethod), pointer :: bcfunc + + ! Register prescribed BC for the solution + bcfunc => hbc1d_Prescribed_burgers1d_shock + call this%hyperbolicBCs%RegisterBoundaryCondition( & + SELF_BC_PRESCRIBED, "prescribed", bcfunc) + + ! Register prescribed BC for the gradient + bcfunc => pbc1d_Prescribed_burgers1d_shock + call this%parabolicBCs%RegisterBoundaryCondition( & + SELF_BC_PRESCRIBED, "prescribed", bcfunc) + endsubroutine + + subroutine hbc1d_Prescribed_burgers1d_shock(bc, mymodel) + class(BoundaryCondition), intent(in) :: bc + class(Model), intent(inout) :: mymodel + integer :: n, iEl, s + real(prec) :: x, jump, spd + + select type (m => mymodel) + class is (burgers1d_shock) + do n = 1, bc%nBoundaries + iEl = bc%elements(n) + s = bc%sides(n) + x = m%geometry%x%boundary(s, iEl, 1) + jump = m%ul - m%ur + spd = 0.5_prec*(m%ul + m%ur) + m%solution%extBoundary(s, iEl, 1) = & + spd - 0.5_prec*tanh((x - spd*m%t - m%x0)*jump & + / (4.0_prec*m%nu)) + enddo + endselect + endsubroutine + + subroutine pbc1d_Prescribed_burgers1d_shock(bc, mymodel) + class(BoundaryCondition), intent(in) :: bc + class(Model), intent(inout) :: mymodel + integer :: n, iEl, s + real(prec) :: x, jump, spd, r, drdx + + select type (m => mymodel) + class is (burgers1d_shock) + do n = 1, bc%nBoundaries + iEl = bc%elements(n) + s = bc%sides(n) + x = m%geometry%x%boundary(s, iEl, 1) + jump = m%ul - m%ur + spd = 0.5_prec*(m%ul + m%ur) + r = (x - spd*m%t - m%x0)*jump/(4.0_prec*m%nu) + drdx = jump/(4.0_prec*m%nu) + m%solutionGradient%extBoundary(s, iEl, 1) = & + -0.5_prec*drdx*(2.0_prec/(exp(r) + exp(-r)))**2 + enddo + endselect + endsubroutine + +endmodule +``` + +### Main Program + +```fortran +program traveling_shock + use self_data + use burgers1d_shock_model + implicit none + + type(burgers1d_shock) :: modelobj + type(Lagrange), target :: interp + type(Mesh1D), target :: mesh + type(Geometry1D), target :: geometry + + ! Create mesh + call mesh%StructuredMesh(nElem=10, x=(/0.0_prec, 1.0_prec/)) + + ! Tag both boundaries as prescribed + call mesh%ResetBoundaryConditionType(SELF_BC_PRESCRIBED, SELF_BC_PRESCRIBED) + + ! Create interpolant and geometry + call interp%Init(N=7, controlNodeType=GAUSS, M=10, targetNodeType=UNIFORM) + call geometry%Init(interp, mesh%nElem) + call geometry%GenerateFromMesh(mesh) + + ! Initialize model (this calls AdditionalInit and MapBoundaryConditions) + call modelobj%Init(mesh, geometry) + modelobj%gradient_enabled = .true. + modelobj%nu = 0.01_prec + + ! Set initial condition, run, clean up... + call modelobj%SetTimeIntegrator('rk3') + call modelobj%ForwardStep(2.0_prec, 1.0e-5_prec, 0.05_prec) + call modelobj%free() + call mesh%free() + call geometry%free() + call interp%free() +endprogram +``` + +## Extending a Built-in Model's BCs + +When your model extends a built-in SELF model that already registers BCs, call the parent's `AdditionalInit` first, then register your additional BCs: + +```fortran +subroutine AdditionalInit_my_extended_model(this) + class(my_extended_model), intent(inout) :: this + procedure(SELF_bcMethod), pointer :: bcfunc + + ! Register parent model's BCs first + call AdditionalInit_LinearEuler2D_t(this) + + ! Then register additional BCs + bcfunc => hbc2d_Prescribed_my_extended_model + call this%hyperbolicBCs%RegisterBoundaryCondition( & + SELF_BC_PRESCRIBED, "prescribed", bcfunc) +endsubroutine +``` + +## Replacing a BC Implementation + +If a BC with a given ID is already registered and you call `RegisterBoundaryCondition` with the same ID again, the procedure pointer is updated to point to your new implementation. The element/side arrays are preserved. This mechanism is used internally by GPU model variants to replace CPU BC routines with GPU-accelerated kernels, but you can also use it to override a parent model's BC behavior. + +## Tips + +- **One BC per ID per list.** Each `bcid` can appear at most once in the hyperbolic list and once in the parabolic list. +- **Mesh tagging must match registration.** If you register a BC for `SELF_BC_PRESCRIBED` but no mesh faces carry that ID, the BC will have `nBoundaries = 0` and its subroutine will never be called. +- **Use `select type` in BC routines.** The `SELF_bcMethod` interface receives `class(Model)`, so you must downcast to access your model's fields. +- **BC routines are called every time step.** Keep them efficient. For GPU models, use device kernels rather than host-side loops. diff --git a/examples/burgers1d_shock.f90 b/examples/burgers1d_shock.f90 index 687300da5..454109008 100644 --- a/examples/burgers1d_shock.f90 +++ b/examples/burgers1d_shock.f90 @@ -39,6 +39,7 @@ module burgers1d_shock_model !! where $s = \frac{u_l + u_r}{2}$ is the shock speed. use self_Burgers1D + use SELF_BoundaryConditions implicit none @@ -49,42 +50,73 @@ module burgers1d_shock_model contains - procedure :: hbc1d_Prescribed => hbc1d_Prescribed_burgers1d_shock ! override for the hyperbolic boundary conditions - procedure :: pbc1d_Prescribed => pbc1d_Prescribed_burgers1d_shock ! override for the parabolic boundary conditions + procedure :: AdditionalInit => AdditionalInit_burgers1d_shock endtype burgers1d_shock contains - pure function hbc1d_Prescribed_burgers1d_shock(this,x,t) result(exts) - class(burgers1d_shock),intent(in) :: this - real(prec),intent(in) :: x - real(prec),intent(in) :: t - real(prec) :: exts(1:this%nvar) + subroutine AdditionalInit_burgers1d_shock(this) + implicit none + class(burgers1d_shock),intent(inout) :: this ! Local - real(prec) :: jump,s - - jump = this%ul-this%ur - s = 0.5_prec*(this%ul+this%ur) - exts(1) = s-0.5_prec*tanh((x-s*t-this%x0)*jump/(4.0_prec*this%nu)) + procedure(SELF_bcMethod),pointer :: bcfunc - endfunction hbc1d_Prescribed_burgers1d_shock + bcfunc => hbc1d_Prescribed_burgers1d_shock + call this%hyperbolicBCs%RegisterBoundaryCondition( & + SELF_BC_PRESCRIBED,"prescribed",bcfunc) - pure function pbc1d_Prescribed_burgers1d_shock(this,x,t) result(extDsdx) - class(burgers1d_shock),intent(in) :: this - real(prec),intent(in) :: x - real(prec),intent(in) :: t - real(prec) :: extDsdx(1:this%nvar) - ! Local - real(prec) :: jump,s,r,drdx + bcfunc => pbc1d_Prescribed_burgers1d_shock + call this%parabolicBCs%RegisterBoundaryCondition( & + SELF_BC_PRESCRIBED,"prescribed",bcfunc) - jump = this%ul-this%ur - s = 0.5_prec*(this%ul+this%ur) - r = (x-s*t-this%x0)*jump/(4.0_prec*this%nu) - drdx = jump/(4.0_prec*this%nu) - extDsdx(1) = -0.5_prec*drdx*(sech(r))**2 + endsubroutine AdditionalInit_burgers1d_shock - endfunction pbc1d_Prescribed_burgers1d_shock + subroutine hbc1d_Prescribed_burgers1d_shock(bc,mymodel) + class(BoundaryCondition),intent(in) :: bc + class(Model),intent(inout) :: mymodel + ! Local + integer :: n,iEl,s + real(prec) :: x,jump,spd + + select type(m => mymodel) + class is(burgers1d_shock) + do n = 1,bc%nBoundaries + iEl = bc%elements(n) + s = bc%sides(n) + x = m%geometry%x%boundary(s,iEl,1) + jump = m%ul-m%ur + spd = 0.5_prec*(m%ul+m%ur) + m%solution%extBoundary(s,iEl,1) = & + spd-0.5_prec*tanh((x-spd*m%t-m%x0)*jump/(4.0_prec*m%nu)) + enddo + endselect + + endsubroutine hbc1d_Prescribed_burgers1d_shock + + subroutine pbc1d_Prescribed_burgers1d_shock(bc,mymodel) + class(BoundaryCondition),intent(in) :: bc + class(Model),intent(inout) :: mymodel + ! Local + integer :: n,iEl,s + real(prec) :: x,jump,spd,r,drdx + + select type(m => mymodel) + class is(burgers1d_shock) + do n = 1,bc%nBoundaries + iEl = bc%elements(n) + s = bc%sides(n) + x = m%geometry%x%boundary(s,iEl,1) + jump = m%ul-m%ur + spd = 0.5_prec*(m%ul+m%ur) + r = (x-spd*m%t-m%x0)*jump/(4.0_prec*m%nu) + drdx = jump/(4.0_prec*m%nu) + m%solutionGradient%extBoundary(s,iEl,1) = & + -0.5_prec*drdx*(sech(r))**2 + enddo + endselect + + endsubroutine pbc1d_Prescribed_burgers1d_shock pure real(prec) function sech(x) result(fx) real(prec),intent(in) :: x diff --git a/examples/linear_euler2d_planewave_propagation.f90 b/examples/linear_euler2d_planewave_propagation.f90 index 85769f197..a412af0f2 100644 --- a/examples/linear_euler2d_planewave_propagation.f90 +++ b/examples/linear_euler2d_planewave_propagation.f90 @@ -36,6 +36,7 @@ module lineareuler2d_planewave_prop_model !! See section 5.4.3 of D. A. Kopriva, "Implementing Spectral Methods for Partial Differential Equations" (2009) use self_lineareuler2d + use SELF_BoundaryConditions implicit none @@ -50,12 +51,28 @@ module lineareuler2d_planewave_prop_model contains procedure :: setInitialCondition - procedure :: hbc2d_Prescribed => hbc2d_Prescribed_lineareuler2d_planewave ! override for the hyperbolic boundary conditions + procedure :: AdditionalInit => AdditionalInit_lineareuler2d_planewave endtype lineareuler2d_planewave contains + subroutine AdditionalInit_lineareuler2d_planewave(this) + implicit none + class(lineareuler2d_planewave),intent(inout) :: this + ! Local + procedure(SELF_bcMethod),pointer :: bcfunc + + ! Register the parent class NoNormalFlow BC first + call AdditionalInit_LinearEuler2D_t(this) + + ! Register the prescribed BC + bcfunc => hbc2d_Prescribed_lineareuler2d_planewave + call this%hyperbolicBCs%RegisterBoundaryCondition( & + SELF_BC_PRESCRIBED,"prescribed",bcfunc) + + endsubroutine AdditionalInit_lineareuler2d_planewave + subroutine setInitialCondition(this) implicit none class(lineareuler2d_planewave),intent(inout) :: this @@ -87,29 +104,39 @@ subroutine setInitialCondition(this) endsubroutine setInitialCondition - pure function hbc2d_Prescribed_lineareuler2d_planewave(this,x,t) result(exts) - class(lineareuler2d_planewave),intent(in) :: this - real(prec),intent(in) :: x(1:2) - real(prec),intent(in) :: t - real(prec) :: exts(1:this%nvar) + subroutine hbc2d_Prescribed_lineareuler2d_planewave(bc,mymodel) + class(BoundaryCondition),intent(in) :: bc + class(Model),intent(inout) :: mymodel ! Local + integer :: n,i,iEl,j + real(prec) :: x(1:2) real(prec) :: p,rho,u,v,phase,shape - p = this%p - rho = this%p/this%c/this%c - u = this%p*this%wx/this%c - v = this%p*this%wy/this%c - - phase = this%wx*(x(1)-this%x0)+this%wy*(x(2)-this%y0)-this%c*this%t - shape = exp(-phase*phase/(this%L*this%L)) - - exts(1) = rho*shape ! density - exts(2) = u*shape ! u - exts(3) = v*shape ! v - exts(4) = p*shape ! pressure - if(.false.) exts(1) = exts(1)+t ! suppress unused-dummy-argument warning - - endfunction hbc2d_Prescribed_lineareuler2d_planewave + select type(m => mymodel) + class is(lineareuler2d_planewave) + p = m%p + rho = m%p/m%c/m%c + u = m%p*m%wx/m%c + v = m%p*m%wy/m%c + + do n = 1,bc%nBoundaries + iEl = bc%elements(n) + j = bc%sides(n) + do i = 1,m%solution%interp%N+1 + x = m%geometry%x%boundary(i,j,iEl,1,1:2) + + phase = m%wx*(x(1)-m%x0)+m%wy*(x(2)-m%y0)-m%c*m%t + shape = exp(-phase*phase/(m%L*m%L)) + + m%solution%extBoundary(i,j,iEl,1) = rho*shape ! density + m%solution%extBoundary(i,j,iEl,2) = u*shape ! u + m%solution%extBoundary(i,j,iEl,3) = v*shape ! v + m%solution%extBoundary(i,j,iEl,4) = p*shape ! pressure + enddo + enddo + endselect + + endsubroutine hbc2d_Prescribed_lineareuler2d_planewave endmodule lineareuler2d_planewave_prop_model diff --git a/examples/linear_euler2d_planewave_reflection.f90 b/examples/linear_euler2d_planewave_reflection.f90 index 508c81c7f..f38f00af9 100644 --- a/examples/linear_euler2d_planewave_reflection.f90 +++ b/examples/linear_euler2d_planewave_reflection.f90 @@ -39,6 +39,7 @@ module lineareuler2d_planewave_model !! on the "eastern" boundary of the domain. use self_lineareuler2d + use SELF_BoundaryConditions implicit none @@ -53,12 +54,28 @@ module lineareuler2d_planewave_model contains procedure :: setInitialCondition - procedure :: hbc2d_Prescribed => hbc2d_Prescribed_lineareuler2d_planewave ! override for the hyperbolic boundary conditions + procedure :: AdditionalInit => AdditionalInit_lineareuler2d_planewave endtype lineareuler2d_planewave contains + subroutine AdditionalInit_lineareuler2d_planewave(this) + implicit none + class(lineareuler2d_planewave),intent(inout) :: this + ! Local + procedure(SELF_bcMethod),pointer :: bcfunc + + ! Register the parent class NoNormalFlow BC first + call AdditionalInit_LinearEuler2D_t(this) + + ! Register the prescribed BC + bcfunc => hbc2d_Prescribed_lineareuler2d_planewave + call this%hyperbolicBCs%RegisterBoundaryCondition( & + SELF_BC_PRESCRIBED,"prescribed",bcfunc) + + endsubroutine AdditionalInit_lineareuler2d_planewave + subroutine setInitialCondition(this) implicit none class(lineareuler2d_planewave),intent(inout) :: this @@ -97,34 +114,44 @@ subroutine setInitialCondition(this) endsubroutine setInitialCondition - pure function hbc2d_Prescribed_lineareuler2d_planewave(this,x,t) result(exts) - class(lineareuler2d_planewave),intent(in) :: this - real(prec),intent(in) :: x(1:2) - real(prec),intent(in) :: t - real(prec) :: exts(1:this%nvar) + subroutine hbc2d_Prescribed_lineareuler2d_planewave(bc,mymodel) + class(BoundaryCondition),intent(in) :: bc + class(Model),intent(inout) :: mymodel ! Local + integer :: n,i,iEl,j + real(prec) :: x(1:2) real(prec) :: p,rho,u,v,phase,shi,shr - p = this%p - rho = this%p/this%c/this%c - u = this%p*this%wx/this%c - v = this%p*this%wy/this%c - - ! Incident wave - phase = this%wx*(x(1)-this%x0)+this%wy*(x(2)-this%y0)-this%c*this%t - shi = exp(-phase*phase/(this%L*this%L)) - - ! Reflected wave - phase = -this%wx*(x(1)+this%x0-2.0_prec)+this%wy*(x(2)-this%y0)-this%c*this%t - shr = exp(-phase*phase/(this%L*this%L)) - - exts(1) = rho*(shi+shr) ! density - exts(2) = u*(shi-shr) ! u - exts(3) = v*(shi+shr) ! v - exts(4) = p*(shi+shr) ! pressure - if(.false.) exts(1) = exts(1)+t ! suppress unused-dummy-argument warning - - endfunction hbc2d_Prescribed_lineareuler2d_planewave + select type(m => mymodel) + class is(lineareuler2d_planewave) + p = m%p + rho = m%p/m%c/m%c + u = m%p*m%wx/m%c + v = m%p*m%wy/m%c + + do n = 1,bc%nBoundaries + iEl = bc%elements(n) + j = bc%sides(n) + do i = 1,m%solution%interp%N+1 + x = m%geometry%x%boundary(i,j,iEl,1,1:2) + + ! Incident wave + phase = m%wx*(x(1)-m%x0)+m%wy*(x(2)-m%y0)-m%c*m%t + shi = exp(-phase*phase/(m%L*m%L)) + + ! Reflected wave + phase = -m%wx*(x(1)+m%x0-2.0_prec)+m%wy*(x(2)-m%y0)-m%c*m%t + shr = exp(-phase*phase/(m%L*m%L)) + + m%solution%extBoundary(i,j,iEl,1) = rho*(shi+shr) ! density + m%solution%extBoundary(i,j,iEl,2) = u*(shi-shr) ! u + m%solution%extBoundary(i,j,iEl,3) = v*(shi+shr) ! v + m%solution%extBoundary(i,j,iEl,4) = p*(shi+shr) ! pressure + enddo + enddo + endselect + + endsubroutine hbc2d_Prescribed_lineareuler2d_planewave endmodule lineareuler2d_planewave_model diff --git a/mkdocs.yml b/mkdocs.yml index b3838f01c..e0d3488fe 100644 --- a/mkdocs.yml +++ b/mkdocs.yml @@ -34,6 +34,7 @@ nav: - Kelvin waves: Tutorials/LinearShallowWater/KelvinWaves.md - Planetary Rossby wave: Tutorials/LinearShallowWater/PlanetaryRossbyWave.md - Models: + - Boundary Conditions: Models/boundary-conditions.md - Viscous Burgers Equation: Models/burgers-equation-model.md - Linear Euler (2D) : Models/linear-euler-2d-model.md - Linear Euler (3D) : Models/linear-euler-3d-model.md @@ -52,6 +53,7 @@ nav: - Provable Stability: Learning/ProvableStability.md - Code: - Software Architecture: Learning/SoftwareArchitecture.md + - Boundary Condition System: Learning/BoundaryConditions.md - Dependencies: Learning/dependencies.md - API Documentation: ford/index.html - Contributing: diff --git a/src/SELF_BoundaryConditions.f90 b/src/SELF_BoundaryConditions.f90 new file mode 100644 index 000000000..59babf4c3 --- /dev/null +++ b/src/SELF_BoundaryConditions.f90 @@ -0,0 +1,215 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +module SELF_BoundaryConditions + + use SELF_SupportRoutines + use SELF_Metadata + use iso_c_binding + + implicit none + integer,parameter :: SELF_BCNAME_LENGTH = 32 + + type BoundaryCondition + procedure(SELF_bcMethod),pointer :: bcMethod => null() + integer :: bcid + character(SELF_BCNAME_LENGTH) :: bcname + integer :: nBoundaries ! Number of boundaries this BC applies to + integer,allocatable :: elements(:) ! List of elements this BC applies to + integer,allocatable :: sides(:) ! List of local sides this BC applies to + type(c_ptr) :: elements_gpu = c_null_ptr ! Device pointer for elements + type(c_ptr) :: sides_gpu = c_null_ptr ! Device pointer for sides + type(BoundaryCondition),pointer :: next => null() + type(BoundaryCondition),pointer :: prev => null() + endtype BoundaryCondition + + type BoundaryConditionList + type(BoundaryCondition),pointer :: current => null() + type(BoundaryCondition),pointer :: head => null() + type(BoundaryCondition),pointer :: tail => null() + integer :: nbc + + contains + procedure,public :: Init => Init_BCList + procedure,public :: Free => Free_BCList + procedure,private :: MoveNext + procedure,private :: rewind + procedure,public :: GetBCForID + generic,public :: RegisterBoundaryCondition => RegisterbcMethod + procedure,private :: RegisterbcMethod + procedure,public :: PopulateBoundaries + + endtype BoundaryConditionList + + interface + subroutine SELF_bcMethod(this,mymodel) + use SELF_Constants,only:prec + use SELF_Model,only:Model + import BoundaryCondition + implicit none + class(BoundaryCondition),intent(in) :: this + class(Model),intent(inout) :: mymodel + endsubroutine SELF_bcMethod + endinterface + +contains + +! //////////////////////////////////////////// ! +! Boundary Condition Methods +! ////////////////////////////////////////////// ! + + subroutine Init_BCList(list) + class(BoundaryConditionList),intent(inout) :: list + list%head => null() + list%tail => null() + list%current => null() + list%nbc = 0 + endsubroutine Init_BCList + + subroutine Free_BCList(list) + class(BoundaryConditionList),intent(inout) :: list + type(BoundaryCondition),pointer :: node,next_node + + node => list%head + do while(associated(node)) + next_node => node%next + nullify(node%bcMethod) + if(allocated(node%elements)) deallocate(node%elements) + if(allocated(node%sides)) deallocate(node%sides) + deallocate(node) + node => next_node + enddo + + call Init_BCList(list) + endsubroutine Free_BCList + + subroutine MoveNext(list) + class(BoundaryConditionList),intent(inout) :: list + if(associated(list%current%next)) then + list%current => list%current%next + else + nullify(list%current) + endif + endsubroutine MoveNext + + subroutine rewind(list) + class(BoundaryConditionList),intent(inout) :: list + list%current => list%head + endsubroutine rewind + + function GetBCForID(list,bcid) result(node) + !! Returns the node associated with the given bcid. + !! If the bcid is not found, a null pointer is returned. + class(BoundaryConditionList),intent(in) :: list + integer,intent(in) :: bcid + type(BoundaryCondition),pointer :: node + + node => list%head + + do while(associated(node)) + if(node%bcid == bcid) then + return + endif + node => node%next + enddo + ! If we reach this point, the bcid was not found + node => null() + + endfunction GetBCForID + + subroutine RegisterbcMethod(list,bcid,bcname,bcfunc) + !! Register a boundary condition function with the given bcid and bcname. + !! If the bcid is already registered, the function pointer is updated. + !! The elements and sides arrays are not allocated here; call + !! PopulateBoundaries after scanning the mesh. + class(BoundaryConditionList),intent(inout) :: list + integer,intent(in) :: bcid + character(*),intent(in) :: bcname + procedure(SELF_bcMethod),pointer,intent(in) :: bcfunc + ! Local + type(BoundaryCondition),pointer :: bc + + ! Check if bcid is registered + bc => list%GetBCForID(bcid) + if(associated(bc)) then + print*,"Boundary condition with ID ",bcid," is already registered." + print*,"Assigning new function to existing BC" + bc%bcMethod => bcfunc + else + allocate(bc) + bc%bcid = bcid + bc%bcname = trim(bcname) + bc%bcMethod => bcfunc + bc%nBoundaries = 0 + nullify(bc%next) + nullify(bc%prev) + + ! Insert at the tail + if(.not. associated(list%head)) then + ! First entry + list%head => bc + list%tail => bc + else + ! Append to tail + bc%prev => list%tail + list%tail%next => bc + list%tail => bc + endif + + list%nbc = list%nbc+1 + list%current => bc + + endif + + endsubroutine RegisterbcMethod + + subroutine PopulateBoundaries(list,bcid,nBoundaries,elements,sides) + !! Populate the elements and sides arrays for a registered boundary condition. + !! Called after scanning the mesh to determine which faces belong to each bcid. + class(BoundaryConditionList),intent(inout) :: list + integer,intent(in) :: bcid + integer,intent(in) :: nBoundaries + integer,intent(in) :: elements(1:nBoundaries) + integer,intent(in) :: sides(1:nBoundaries) + ! Local + type(BoundaryCondition),pointer :: bc + + bc => list%GetBCForID(bcid) + if(.not. associated(bc)) then + return + endif + + bc%nBoundaries = nBoundaries + if(allocated(bc%elements)) deallocate(bc%elements) + if(allocated(bc%sides)) deallocate(bc%sides) + allocate(bc%elements(1:nBoundaries)) + allocate(bc%sides(1:nBoundaries)) + bc%elements = elements + bc%sides = sides + + endsubroutine PopulateBoundaries + +endmodule SELF_BoundaryConditions diff --git a/src/SELF_DGModel1D_t.f90 b/src/SELF_DGModel1D_t.f90 index 75122257d..d7937c260 100644 --- a/src/SELF_DGModel1D_t.f90 +++ b/src/SELF_DGModel1D_t.f90 @@ -34,6 +34,7 @@ module SELF_DGModel1D_t use HDF5 use FEQParse use SELF_Model + use SELF_BoundaryConditions implicit none @@ -47,12 +48,15 @@ module SELF_DGModel1D_t type(MappedScalar1D) :: workSol type(Mesh1D),pointer :: mesh type(Geometry1D),pointer :: geometry + type(BoundaryConditionList) :: hyperbolicBCs + type(BoundaryConditionList) :: parabolicBCs contains procedure :: Init => Init_DGModel1D_t procedure :: SetMetadata => SetMetadata_DGModel1D_t procedure :: Free => Free_DGModel1D_t + procedure :: MapBoundaryConditions => MapBoundaryConditions_DGModel1D_t procedure :: CalculateEntropy => CalculateEntropy_DGModel1D_t procedure :: BoundaryFlux => BoundaryFlux_DGModel1D_t @@ -106,8 +110,13 @@ subroutine Init_DGModel1D_t(this,mesh,geometry) call this%flux%AssociateGeometry(geometry) call this%fluxDivergence%AssociateGeometry(geometry) + call this%hyperbolicBCs%Init() + call this%parabolicBCs%Init() + call this%AdditionalInit() + call this%MapBoundaryConditions() + call this%SetMetadata() endsubroutine Init_DGModel1D_t @@ -145,6 +154,8 @@ subroutine Free_DGModel1D_t(this) call this%flux%Free() call this%source%Free() call this%fluxDivergence%Free() + call this%hyperbolicBCs%Free() + call this%parabolicBCs%Free() call this%AdditionalFree() endsubroutine Free_DGModel1D_t @@ -315,127 +326,125 @@ subroutine CalculateEntropy_DGModel1D_t(this) endsubroutine CalculateEntropy_DGModel1D_t - subroutine setboundarycondition_DGModel1D_t(this) - ! Here, we use the pre-tendency method to calculate the - ! derivative of the solution using a bassi-rebay method - ! We then do a boundary interpolation and side exchange - ! on the gradient field + subroutine MapBoundaryConditions_DGModel1D_t(this) + !! Scan the mesh boundary condition IDs and populate the elements/sides + !! arrays for each registered boundary condition. implicit none class(DGModel1D_t),intent(inout) :: this - ! local - integer :: N,nelem - real(prec) :: x - - nelem = this%geometry%nelem ! number of elements in the mesh - N = this%solution%interp%N ! polynomial degree - - ! left-most boundary - if(this%mesh%bcid(1) == SELF_BC_PRESCRIBED) then - - x = this%geometry%x%boundary(1,1,1) - this%solution%extBoundary(1,1,1:this%nvar) = & - this%hbc1d_Prescribed(x,this%t) - - elseif(this%mesh%bcid(1) == SELF_BC_RADIATION) then - - this%solution%extBoundary(1,1,1:this%nvar) = & - this%hbc1d_Radiation(this%solution%boundary(1,1,1:this%nvar),-1.0_prec) - - elseif(this%mesh%bcid(1) == SELF_BC_NONORMALFLOW) then - - this%solution%extBoundary(1,1,1:this%nvar) = & - this%hbc1d_NoNormalFlow(this%solution%boundary(1,1,1:this%nvar),-1.0_prec) - - else ! Periodic - - this%solution%extBoundary(1,1,1:this%nvar) = this%solution%boundary(2,nelem,1:this%nvar) - - endif - - ! right-most boundary - if(this%mesh%bcid(1) == SELF_BC_PRESCRIBED) then - - x = this%geometry%x%boundary(2,nelem,1) - this%solution%extBoundary(2,nelem,1:this%nvar) = & - this%hbc1d_Prescribed(x,this%t) - - elseif(this%mesh%bcid(1) == SELF_BC_RADIATION) then + ! Local + type(BoundaryCondition),pointer :: bc + integer :: nelem,count,n + integer :: elems(2),sds(2) + + nelem = this%mesh%nElem + + ! Map hyperbolic BCs + bc => this%hyperbolicBCs%head + do while(associated(bc)) + count = 0 + if(this%mesh%bcid(1) == bc%bcid) count = count+1 + if(this%mesh%bcid(2) == bc%bcid) count = count+1 + + if(count > 0) then + n = 0 + if(this%mesh%bcid(1) == bc%bcid) then + n = n+1 + elems(n) = 1 + sds(n) = 1 + endif + if(this%mesh%bcid(2) == bc%bcid) then + n = n+1 + elems(n) = nelem + sds(n) = 2 + endif + call this%hyperbolicBCs%PopulateBoundaries(bc%bcid,count, & + elems(1:count),sds(1:count)) + endif + bc => bc%next + enddo - this%solution%extBoundary(2,nelem,1:this%nvar) = & - this%hbc1d_Radiation(this%solution%boundary(2,nelem,1:this%nvar),-1.0_prec) + ! Map parabolic BCs + bc => this%parabolicBCs%head + do while(associated(bc)) + count = 0 + if(this%mesh%bcid(1) == bc%bcid) count = count+1 + if(this%mesh%bcid(2) == bc%bcid) count = count+1 + + if(count > 0) then + n = 0 + if(this%mesh%bcid(1) == bc%bcid) then + n = n+1 + elems(n) = 1 + sds(n) = 1 + endif + if(this%mesh%bcid(2) == bc%bcid) then + n = n+1 + elems(n) = nelem + sds(n) = 2 + endif + call this%parabolicBCs%PopulateBoundaries(bc%bcid,count, & + elems(1:count),sds(1:count)) + endif + bc => bc%next + enddo - elseif(this%mesh%bcid(1) == SELF_BC_NONORMALFLOW) then + endsubroutine MapBoundaryConditions_DGModel1D_t - this%solution%extBoundary(2,nelem,1:this%nvar) = & - this%hbc1d_NoNormalFlow(this%solution%boundary(2,nelem,1:this%nvar),-1.0_prec) + subroutine setboundarycondition_DGModel1D_t(this) + !! Apply boundary conditions for the solution. + !! Periodic boundaries are set as the default; registered + !! boundary conditions overwrite specific endpoints. + implicit none + class(DGModel1D_t),intent(inout) :: this + ! Local + type(BoundaryCondition),pointer :: bc + procedure(SELF_bcMethod),pointer :: apply_bc + integer :: nelem - else ! Periodic + nelem = this%mesh%nElem - this%solution%extBoundary(2,nelem,1:this%nvar) = this%solution%boundary(1,1,1:this%nvar) + ! Default: periodic boundary conditions + this%solution%extBoundary(1,1,1:this%nvar) = & + this%solution%boundary(2,nelem,1:this%nvar) + this%solution%extBoundary(2,nelem,1:this%nvar) = & + this%solution%boundary(1,1,1:this%nvar) - endif + ! Apply registered boundary conditions + bc => this%hyperbolicBCs%head + do while(associated(bc)) + apply_bc => bc%bcMethod + call apply_bc(bc,this) + bc => bc%next + enddo endsubroutine setboundarycondition_DGModel1D_t subroutine setgradientboundarycondition_DGModel1D_t(this) - ! Here, we set the boundary conditions for the - ! solution and the solution gradient at the left - ! and right most boundaries. - ! - ! Here, we use periodic boundary conditions + !! Apply boundary conditions for the solution gradient. + !! Periodic boundaries are set as the default; registered + !! boundary conditions overwrite specific endpoints. implicit none class(DGModel1D_t),intent(inout) :: this - ! local - real(prec) :: x + ! Local + type(BoundaryCondition),pointer :: bc + procedure(SELF_bcMethod),pointer :: apply_bc integer :: nelem - nelem = this%geometry%nelem ! number of elements in the mesh - - ! left-most boundary - if(this%mesh%bcid(1) == SELF_BC_PRESCRIBED) then - - x = this%geometry%x%boundary(1,1,1) - this%solutionGradient%extBoundary(1,1,1:this%nvar) = & - this%pbc1d_Prescribed(x,this%t) - - elseif(this%mesh%bcid(1) == SELF_BC_RADIATION) then - - this%solutionGradient%extBoundary(1,1,1:this%nvar) = & - this%pbc1d_Radiation(this%solutionGradient%boundary(1,1,1:this%nvar),-1.0_prec) - - elseif(this%mesh%bcid(1) == SELF_BC_NONORMALFLOW) then - - this%solutionGradient%extBoundary(1,1,1:this%nvar) = & - this%pbc1d_NoNormalFlow(this%solutionGradient%boundary(1,1,1:this%nvar),-1.0_prec) - - else ! Periodic + nelem = this%mesh%nElem - this%solutionGradient%extBoundary(1,1,1:this%nvar) = this%solutionGradient%boundary(2,nelem,1:this%nvar) + ! Default: periodic boundary conditions + this%solutionGradient%extBoundary(1,1,1:this%nvar) = & + this%solutionGradient%boundary(2,nelem,1:this%nvar) + this%solutionGradient%extBoundary(2,nelem,1:this%nvar) = & + this%solutionGradient%boundary(1,1,1:this%nvar) - endif - - ! right-most boundary - if(this%mesh%bcid(1) == SELF_BC_PRESCRIBED) then - - x = this%geometry%x%boundary(2,nelem,1) - this%solutionGradient%extBoundary(2,nelem,1:this%nvar) = & - this%pbc1d_Prescribed(x,this%t) - - elseif(this%mesh%bcid(1) == SELF_BC_RADIATION) then - - this%solutionGradient%extBoundary(2,nelem,1:this%nvar) = & - this%pbc1d_Radiation(this%solutionGradient%boundary(2,nelem,1:this%nvar),-1.0_prec) - - elseif(this%mesh%bcid(1) == SELF_BC_NONORMALFLOW) then - - this%solutionGradient%extBoundary(2,nelem,1:this%nvar) = & - this%pbc1d_NoNormalFlow(this%solutionGradient%boundary(2,nelem,1:this%nvar),-1.0_prec) - - else ! Periodic - - this%solutionGradient%extBoundary(2,nelem,1:this%nvar) = this%solutionGradient%boundary(1,1,1:this%nvar) - - endif + ! Apply registered boundary conditions + bc => this%parabolicBCs%head + do while(associated(bc)) + apply_bc => bc%bcMethod + call apply_bc(bc,this) + bc => bc%next + enddo endsubroutine setgradientboundarycondition_DGModel1D_t diff --git a/src/SELF_DGModel2D_t.f90 b/src/SELF_DGModel2D_t.f90 index d87b19a05..01a237dda 100644 --- a/src/SELF_DGModel2D_t.f90 +++ b/src/SELF_DGModel2D_t.f90 @@ -36,6 +36,7 @@ module SELF_DGModel2D_t use HDF5 use FEQParse use SELF_Model + use SELF_BoundaryConditions implicit none @@ -49,12 +50,15 @@ module SELF_DGModel2D_t type(MappedScalar2D) :: workSol type(Mesh2D),pointer :: mesh type(SEMQuad),pointer :: geometry + type(BoundaryConditionList) :: hyperbolicBCs + type(BoundaryConditionList) :: parabolicBCs contains procedure :: Init => Init_DGModel2D_t procedure :: SetMetadata => SetMetadata_DGModel2D_t procedure :: Free => Free_DGModel2D_t + procedure :: MapBoundaryConditions => MapBoundaryConditions_DGModel2D_t procedure :: CalculateEntropy => CalculateEntropy_DGModel2D_t procedure :: BoundaryFlux => BoundaryFlux_DGModel2D_t @@ -109,8 +113,13 @@ subroutine Init_DGModel2D_t(this,mesh,geometry) call this%flux%AssociateGeometry(geometry) call this%fluxDivergence%AssociateGeometry(geometry) + call this%hyperbolicBCs%Init() + call this%parabolicBCs%Init() + call this%AdditionalInit() + call this%MapBoundaryConditions() + call this%SetMetadata() endsubroutine Init_DGModel2D_t @@ -143,6 +152,8 @@ subroutine Free_DGModel2D_t(this) call this%flux%Free() call this%source%Free() call this%fluxDivergence%Free() + call this%hyperbolicBCs%Free() + call this%parabolicBCs%Free() call this%AdditionalFree() endsubroutine Free_DGModel2D_t @@ -434,107 +445,119 @@ subroutine sourcemethod_DGModel2D_t(this) endsubroutine sourcemethod_DGModel2D_t - subroutine setboundarycondition_DGModel2D_t(this) - !! Boundary conditions for the solution are set to - !! 0 for the external state to provide radiation type - !! boundary conditions. + subroutine MapBoundaryConditions_DGModel2D_t(this) + !! Scan the mesh sideInfo and populate the elements/sides + !! arrays for each registered boundary condition. implicit none class(DGModel2D_t),intent(inout) :: this - ! local - integer :: i,iEl,j,e2,bcid - real(prec) :: nhat(1:2),x(1:2) - - do concurrent(j=1:4,iel=1:this%mesh%nElem) - - bcid = this%mesh%sideInfo(5,j,iEl) ! Boundary Condition ID - e2 = this%mesh%sideInfo(3,j,iEl) ! Neighboring Element ID - - if(e2 == 0) then - if(bcid == SELF_BC_PRESCRIBED) then - - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - x = this%geometry%x%boundary(i,j,iEl,1,1:2) - - this%solution%extBoundary(i,j,iEl,1:this%nvar) = & - this%hbc2d_Prescribed(x,this%t) - enddo - - elseif(bcid == SELF_BC_RADIATION) then - - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) + ! Local + type(BoundaryCondition),pointer :: bc + integer :: iEl,j,e2,bcid + integer :: count,n + integer,allocatable :: elems(:),sds(:) + + ! Map hyperbolic BCs + bc => this%hyperbolicBCs%head + do while(associated(bc)) + ! Pass 1: count boundary faces for this bcid + count = 0 + do iEl = 1,this%mesh%nElem + do j = 1,4 + e2 = this%mesh%sideInfo(3,j,iEl) + bcid = this%mesh%sideInfo(5,j,iEl) + if(e2 == 0 .and. bcid == bc%bcid) count = count+1 + enddo + enddo - this%solution%extBoundary(i,j,iEl,1:this%nvar) = & - this%hbc2d_Radiation(this%solution%boundary(i,j,iEl,1:this%nvar),nhat) + if(count > 0) then + ! Pass 2: fill element/side arrays + allocate(elems(count),sds(count)) + n = 0 + do iEl = 1,this%mesh%nElem + do j = 1,4 + e2 = this%mesh%sideInfo(3,j,iEl) + bcid = this%mesh%sideInfo(5,j,iEl) + if(e2 == 0 .and. bcid == bc%bcid) then + n = n+1 + elems(n) = iEl + sds(n) = j + endif enddo + enddo + call this%hyperbolicBCs%PopulateBoundaries(bc%bcid,count,elems,sds) + deallocate(elems,sds) + endif + bc => bc%next + enddo - elseif(bcid == SELF_BC_NONORMALFLOW) then - - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) + ! Map parabolic BCs + bc => this%parabolicBCs%head + do while(associated(bc)) + count = 0 + do iEl = 1,this%mesh%nElem + do j = 1,4 + e2 = this%mesh%sideInfo(3,j,iEl) + bcid = this%mesh%sideInfo(5,j,iEl) + if(e2 == 0 .and. bcid == bc%bcid) count = count+1 + enddo + enddo - this%solution%extBoundary(i,j,iEl,1:this%nvar) = & - this%hbc2d_NoNormalFlow(this%solution%boundary(i,j,iEl,1:this%nvar),nhat) + if(count > 0) then + allocate(elems(count),sds(count)) + n = 0 + do iEl = 1,this%mesh%nElem + do j = 1,4 + e2 = this%mesh%sideInfo(3,j,iEl) + bcid = this%mesh%sideInfo(5,j,iEl) + if(e2 == 0 .and. bcid == bc%bcid) then + n = n+1 + elems(n) = iEl + sds(n) = j + endif enddo - - endif + enddo + call this%parabolicBCs%PopulateBoundaries(bc%bcid,count,elems,sds) + deallocate(elems,sds) endif + bc => bc%next + enddo + + endsubroutine MapBoundaryConditions_DGModel2D_t + subroutine setboundarycondition_DGModel2D_t(this) + !! Apply registered boundary conditions for the solution. + !! Each boundary condition method loops over its own + !! boundary faces. + implicit none + class(DGModel2D_t),intent(inout) :: this + ! Local + type(BoundaryCondition),pointer :: bc + procedure(SELF_bcMethod),pointer :: apply_bc + + bc => this%hyperbolicBCs%head + do while(associated(bc)) + apply_bc => bc%bcMethod + call apply_bc(bc,this) + bc => bc%next enddo endsubroutine setboundarycondition_DGModel2D_t subroutine setgradientboundarycondition_DGModel2D_t(this) - !! Boundary conditions for the solution are set to - !! 0 for the external state to provide radiation type - !! boundary conditions. + !! Apply registered boundary conditions for the solution gradient. + !! Each boundary condition method loops over its own + !! boundary faces. implicit none class(DGModel2D_t),intent(inout) :: this - ! local - integer :: i,iEl,j,e2,bcid - real(prec) :: dsdx(1:this%nvar,1:2) - real(prec) :: nhat(1:2),x(1:2) - - do concurrent(j=1:4,iel=1:this%mesh%nElem) - - bcid = this%mesh%sideInfo(5,j,iEl) ! Boundary Condition ID - e2 = this%mesh%sideInfo(3,j,iEl) ! Neighboring Element ID - - if(e2 == 0) then - if(bcid == SELF_BC_PRESCRIBED) then - - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - x = this%geometry%x%boundary(i,j,iEl,1,1:2) - - this%solutiongradient%extBoundary(i,j,iEl,1:this%nvar,1:2) = & - this%pbc2d_Prescribed(x,this%t) - enddo - - elseif(bcid == SELF_BC_RADIATION) then - - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) - - dsdx = this%solutiongradient%boundary(i,j,iEl,1:this%nvar,1:2) - - this%solutiongradient%extBoundary(i,j,iEl,1:this%nvar,1:2) = & - this%pbc2d_Radiation(dsdx,nhat) - enddo - - elseif(bcid == SELF_BC_NONORMALFLOW) then - - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) - - dsdx = this%solutiongradient%boundary(i,j,iEl,1:this%nvar,1:2) - - this%solutiongradient%extBoundary(i,j,iEl,1:this%nvar,1:2) = & - this%pbc2d_NoNormalFlow(dsdx,nhat) - enddo - - endif - endif - + ! Local + type(BoundaryCondition),pointer :: bc + procedure(SELF_bcMethod),pointer :: apply_bc + + bc => this%parabolicBCs%head + do while(associated(bc)) + apply_bc => bc%bcMethod + call apply_bc(bc,this) + bc => bc%next enddo endsubroutine setgradientboundarycondition_DGModel2D_t diff --git a/src/SELF_DGModel3D_t.f90 b/src/SELF_DGModel3D_t.f90 index f0d9c3de1..e376937e4 100644 --- a/src/SELF_DGModel3D_t.f90 +++ b/src/SELF_DGModel3D_t.f90 @@ -36,6 +36,7 @@ module SELF_DGModel3D_t use HDF5 use FEQParse use SELF_Model + use SELF_BoundaryConditions implicit none @@ -49,12 +50,15 @@ module SELF_DGModel3D_t type(MappedScalar3D) :: workSol type(Mesh3D),pointer :: mesh type(SEMHex),pointer :: geometry + type(BoundaryConditionList) :: hyperbolicBCs + type(BoundaryConditionList) :: parabolicBCs contains procedure :: Init => Init_DGModel3D_t procedure :: SetMetadata => SetMetadata_DGModel3D_t procedure :: Free => Free_DGModel3D_t + procedure :: MapBoundaryConditions => MapBoundaryConditions_DGModel3D_t procedure :: CalculateEntropy => CalculateEntropy_DGModel3D_t procedure :: BoundaryFlux => BoundaryFlux_DGModel3D_t @@ -109,8 +113,13 @@ subroutine Init_DGModel3D_t(this,mesh,geometry) call this%flux%AssociateGeometry(geometry) call this%fluxDivergence%AssociateGeometry(geometry) + call this%hyperbolicBCs%Init() + call this%parabolicBCs%Init() + call this%AdditionalInit() + call this%MapBoundaryConditions() + call this%SetMetadata() endsubroutine Init_DGModel3D_t @@ -143,6 +152,8 @@ subroutine Free_DGModel3D_t(this) call this%flux%Free() call this%source%Free() call this%fluxDivergence%Free() + call this%hyperbolicBCs%Free() + call this%parabolicBCs%Free() call this%AdditionalFree() endsubroutine Free_DGModel3D_t @@ -431,119 +442,113 @@ subroutine sourcemethod_DGModel3D_t(this) endsubroutine sourcemethod_DGModel3D_t - subroutine setboundarycondition_DGModel3D_t(this) - !! Boundary conditions for the solution are set to - !! 0 for the external state to provide radiation type - !! boundary conditions. + subroutine MapBoundaryConditions_DGModel3D_t(this) + !! Scan the mesh sideInfo and populate the elements/sides + !! arrays for each registered boundary condition. implicit none class(DGModel3D_t),intent(inout) :: this - ! local - integer :: i,iEl,j,k,e2,bcid - real(prec) :: nhat(1:3),x(1:3) - - do concurrent(k=1:6,iel=1:this%mesh%nElem) - - bcid = this%mesh%sideInfo(5,k,iEl) ! Boundary Condition ID - e2 = this%mesh%sideInfo(3,k,iEl) ! Neighboring Element ID - - if(e2 == 0) then - if(bcid == SELF_BC_PRESCRIBED) then - - do j = 1,this%solution%interp%N+1 ! Loop over quadrature points - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - x = this%geometry%x%boundary(i,j,k,iEl,1,1:3) - - this%solution%extBoundary(i,j,k,iEl,1:this%nvar) = & - this%hbc3d_Prescribed(x,this%t) - enddo - enddo - - elseif(bcid == SELF_BC_RADIATION) then - - do j = 1,this%solution%interp%N+1 ! Loop over quadrature points - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) + ! Local + type(BoundaryCondition),pointer :: bc + integer :: iEl,k,e2,bcid + integer :: count,n + integer,allocatable :: elems(:),sds(:) + + ! Map hyperbolic BCs + bc => this%hyperbolicBCs%head + do while(associated(bc)) + count = 0 + do iEl = 1,this%mesh%nElem + do k = 1,6 + e2 = this%mesh%sideInfo(3,k,iEl) + bcid = this%mesh%sideInfo(5,k,iEl) + if(e2 == 0 .and. bcid == bc%bcid) count = count+1 + enddo + enddo - this%solution%extBoundary(i,j,k,iEl,1:this%nvar) = & - this%hbc3d_Radiation(this%solution%boundary(i,j,k,iEl,1:this%nvar),nhat) - enddo + if(count > 0) then + allocate(elems(count),sds(count)) + n = 0 + do iEl = 1,this%mesh%nElem + do k = 1,6 + e2 = this%mesh%sideInfo(3,k,iEl) + bcid = this%mesh%sideInfo(5,k,iEl) + if(e2 == 0 .and. bcid == bc%bcid) then + n = n+1 + elems(n) = iEl + sds(n) = k + endif enddo + enddo + call this%hyperbolicBCs%PopulateBoundaries(bc%bcid,count,elems,sds) + deallocate(elems,sds) + endif + bc => bc%next + enddo - elseif(bcid == SELF_BC_NONORMALFLOW) then - - do j = 1,this%solution%interp%N+1 ! Loop over quadrature points - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) + ! Map parabolic BCs + bc => this%parabolicBCs%head + do while(associated(bc)) + count = 0 + do iEl = 1,this%mesh%nElem + do k = 1,6 + e2 = this%mesh%sideInfo(3,k,iEl) + bcid = this%mesh%sideInfo(5,k,iEl) + if(e2 == 0 .and. bcid == bc%bcid) count = count+1 + enddo + enddo - this%solution%extBoundary(i,j,k,iEl,1:this%nvar) = & - this%hbc3d_NoNormalFlow(this%solution%boundary(i,j,k,iEl,1:this%nvar),nhat) - enddo + if(count > 0) then + allocate(elems(count),sds(count)) + n = 0 + do iEl = 1,this%mesh%nElem + do k = 1,6 + e2 = this%mesh%sideInfo(3,k,iEl) + bcid = this%mesh%sideInfo(5,k,iEl) + if(e2 == 0 .and. bcid == bc%bcid) then + n = n+1 + elems(n) = iEl + sds(n) = k + endif enddo - - endif + enddo + call this%parabolicBCs%PopulateBoundaries(bc%bcid,count,elems,sds) + deallocate(elems,sds) endif + bc => bc%next + enddo + + endsubroutine MapBoundaryConditions_DGModel3D_t + subroutine setboundarycondition_DGModel3D_t(this) + !! Apply registered boundary conditions for the solution. + implicit none + class(DGModel3D_t),intent(inout) :: this + ! Local + type(BoundaryCondition),pointer :: bc + procedure(SELF_bcMethod),pointer :: apply_bc + + bc => this%hyperbolicBCs%head + do while(associated(bc)) + apply_bc => bc%bcMethod + call apply_bc(bc,this) + bc => bc%next enddo endsubroutine setboundarycondition_DGModel3D_t subroutine setgradientboundarycondition_DGModel3D_t(this) - !! Boundary conditions for the solution are set to - !! 0 for the external state to provide radiation type - !! boundary conditions. + !! Apply registered boundary conditions for the solution gradient. implicit none class(DGModel3D_t),intent(inout) :: this - ! local - integer :: i,iEl,j,k,e2,bcid - real(prec) :: dsdx(1:this%nvar,1:3) - real(prec) :: nhat(1:3),x(1:3) - - do concurrent(k=1:6,iel=1:this%mesh%nElem) - - bcid = this%mesh%sideInfo(5,k,iEl) ! Boundary Condition ID - e2 = this%mesh%sideInfo(3,k,iEl) ! Neighboring Element ID - - if(e2 == 0) then - if(bcid == SELF_BC_PRESCRIBED) then - - do j = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - x = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) - - this%solutiongradient%extBoundary(i,j,k,iEl,1:this%nvar,1:3) = & - this%pbc3d_Prescribed(x,this%t) - enddo - enddo - - elseif(bcid == SELF_BC_RADIATION) then - - do j = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) - - dsdx = this%solutiongradient%boundary(i,j,k,iEl,1:this%nvar,1:3) - - this%solutiongradient%extBoundary(i,j,k,iEl,1:this%nvar,1:3) = & - this%pbc3d_Radiation(dsdx,nhat) - enddo - enddo - - elseif(bcid == SELF_BC_NONORMALFLOW) then - - do j = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) - - dsdx = this%solutiongradient%boundary(i,j,k,iEl,1:this%nvar,1:3) - - this%solutiongradient%extBoundary(i,j,k,iEl,1:this%nvar,1:3) = & - this%pbc3d_NoNormalFlow(dsdx,nhat) - enddo - enddo - - endif - endif - + ! Local + type(BoundaryCondition),pointer :: bc + procedure(SELF_bcMethod),pointer :: apply_bc + + bc => this%parabolicBCs%head + do while(associated(bc)) + apply_bc => bc%bcMethod + call apply_bc(bc,this) + bc => bc%next enddo endsubroutine setgradientboundarycondition_DGModel3D_t diff --git a/src/SELF_ECAdvection2D_t.f90 b/src/SELF_ECAdvection2D_t.f90 index 807df153f..e59ae6008 100644 --- a/src/SELF_ECAdvection2D_t.f90 +++ b/src/SELF_ECAdvection2D_t.f90 @@ -43,6 +43,7 @@ module SELF_ECAdvection2D_t use SELF_ECDGModel2D_t use SELF_mesh + use SELF_BoundaryConditions implicit none @@ -56,11 +57,7 @@ module SELF_ECAdvection2D_t procedure :: entropy_func => entropy_func_ECAdvection2D_t procedure :: twopointflux2d => twopointflux2d_ECAdvection2D_t procedure :: riemannflux2d => riemannflux2d_ECAdvection2D_t - - ! Mirror BC: sets extBoundary = interior boundary state so that - ! sR = sL at all domain faces. Used to eliminate boundary dissipation - ! when testing entropy conservation. - procedure :: hbc2d_NoNormalFlow => hbc2d_NoNormalFlow_ECAdvection2D_t + procedure :: AdditionalInit => AdditionalInit_ECAdvection2D_t endtype ECAdvection2D_t @@ -113,19 +110,40 @@ pure function riemannflux2d_ECAdvection2D_t(this,sL,sR,dsdx,nhat) result(flux) endfunction riemannflux2d_ECAdvection2D_t - pure function hbc2d_NoNormalFlow_ECAdvection2D_t(this,s,nhat) result(exts) + subroutine AdditionalInit_ECAdvection2D_t(this) + implicit none + class(ECAdvection2D_t),intent(inout) :: this + ! Local + procedure(SELF_bcMethod),pointer :: bcfunc + + bcfunc => hbc2d_NoNormalFlow_ECAdvection2D + call this%hyperbolicBCs%RegisterBoundaryCondition( & + SELF_BC_NONORMALFLOW,"no_normal_flow",bcfunc) + + endsubroutine AdditionalInit_ECAdvection2D_t + + subroutine hbc2d_NoNormalFlow_ECAdvection2D(bc,mymodel) !! Mirror boundary condition: sets extBoundary = interior state. !! With the LLF Riemann flux, this gives sR = sL at the boundary, !! so the Riemann flux reduces to the central flux (a.n)*s — no !! dissipation. Use this BC when testing entropy conservation. - class(ECAdvection2D_t),intent(in) :: this - real(prec),intent(in) :: s(1:this%nvar) - real(prec),intent(in) :: nhat(1:2) - real(prec) :: exts(1:this%nvar) - - exts = s - if(.false.) exts(1) = exts(1)+nhat(1) ! suppress unused-dummy-argument warning - - endfunction hbc2d_NoNormalFlow_ECAdvection2D_t + class(BoundaryCondition),intent(in) :: bc + class(Model),intent(inout) :: mymodel + ! Local + integer :: n,i,iEl,j + + select type(m => mymodel) + class is(ECAdvection2D_t) + do n = 1,bc%nBoundaries + iEl = bc%elements(n) + j = bc%sides(n) + do i = 1,m%solution%interp%N+1 + m%solution%extBoundary(i,j,iEl,1:m%nvar) = & + m%solution%boundary(i,j,iEl,1:m%nvar) + enddo + enddo + endselect + + endsubroutine hbc2d_NoNormalFlow_ECAdvection2D endmodule SELF_ECAdvection2D_t diff --git a/src/SELF_ECAdvection3D_t.f90 b/src/SELF_ECAdvection3D_t.f90 index 4247992c5..96f787c0b 100644 --- a/src/SELF_ECAdvection3D_t.f90 +++ b/src/SELF_ECAdvection3D_t.f90 @@ -41,6 +41,7 @@ module SELF_ECAdvection3D_t use SELF_ECDGModel3D_t use SELF_mesh + use SELF_BoundaryConditions implicit none @@ -55,7 +56,7 @@ module SELF_ECAdvection3D_t procedure :: entropy_func => entropy_func_ECAdvection3D_t procedure :: twopointflux3d => twopointflux3d_ECAdvection3D_t procedure :: riemannflux3d => riemannflux3d_ECAdvection3D_t - procedure :: hbc3d_NoNormalFlow => hbc3d_NoNormalFlow_ECAdvection3D_t + procedure :: AdditionalInit => AdditionalInit_ECAdvection3D_t endtype ECAdvection3D_t @@ -108,16 +109,39 @@ pure function riemannflux3d_ECAdvection3D_t(this,sL,sR,dsdx,nhat) result(flux) endfunction riemannflux3d_ECAdvection3D_t - pure function hbc3d_NoNormalFlow_ECAdvection3D_t(this,s,nhat) result(exts) - !! Mirror boundary condition: sets extBoundary = interior state. - class(ECAdvection3D_t),intent(in) :: this - real(prec),intent(in) :: s(1:this%nvar) - real(prec),intent(in) :: nhat(1:3) - real(prec) :: exts(1:this%nvar) + subroutine AdditionalInit_ECAdvection3D_t(this) + implicit none + class(ECAdvection3D_t),intent(inout) :: this + ! Local + procedure(SELF_bcMethod),pointer :: bcfunc + + bcfunc => hbc3d_NoNormalFlow_ECAdvection3D + call this%hyperbolicBCs%RegisterBoundaryCondition( & + SELF_BC_NONORMALFLOW,"no_normal_flow",bcfunc) - exts = s - if(.false.) exts(1) = exts(1)+nhat(1) ! suppress unused-dummy-argument warning + endsubroutine AdditionalInit_ECAdvection3D_t - endfunction hbc3d_NoNormalFlow_ECAdvection3D_t + subroutine hbc3d_NoNormalFlow_ECAdvection3D(bc,mymodel) + !! Mirror boundary condition: sets extBoundary = interior state. + class(BoundaryCondition),intent(in) :: bc + class(Model),intent(inout) :: mymodel + ! Local + integer :: n,i,j,iEl,k + + select type(m => mymodel) + class is(ECAdvection3D_t) + do n = 1,bc%nBoundaries + iEl = bc%elements(n) + k = bc%sides(n) + do j = 1,m%solution%interp%N+1 + do i = 1,m%solution%interp%N+1 + m%solution%extBoundary(i,j,k,iEl,1:m%nvar) = & + m%solution%boundary(i,j,k,iEl,1:m%nvar) + enddo + enddo + enddo + endselect + + endsubroutine hbc3d_NoNormalFlow_ECAdvection3D endmodule SELF_ECAdvection3D_t diff --git a/src/SELF_LinearEuler2D_t.f90 b/src/SELF_LinearEuler2D_t.f90 index b59ad51a2..ba757af09 100644 --- a/src/SELF_LinearEuler2D_t.f90 +++ b/src/SELF_LinearEuler2D_t.f90 @@ -57,6 +57,7 @@ module self_LinearEuler2D_t use self_model use self_dgmodel2d use self_mesh + use SELF_BoundaryConditions implicit none @@ -69,8 +70,8 @@ module self_LinearEuler2D_t contains procedure :: SetNumberOfVariables => SetNumberOfVariables_LinearEuler2D_t procedure :: SetMetadata => SetMetadata_LinearEuler2D_t + procedure :: AdditionalInit => AdditionalInit_LinearEuler2D_t procedure :: entropy_func => entropy_func_LinearEuler2D_t - procedure :: hbc2d_NoNormalFlow => hbc2d_NoNormalFlow_LinearEuler2D_t procedure :: flux2d => flux2d_LinearEuler2D_t procedure :: riemannflux2d => riemannflux2d_LinearEuler2D_t !procedure :: source2d => source2d_LinearEuler2D_t @@ -121,18 +122,47 @@ pure function entropy_func_LinearEuler2D_t(this,s) result(e) endfunction entropy_func_LinearEuler2D_t - pure function hbc2d_NoNormalFlow_LinearEuler2D_t(this,s,nhat) result(exts) - class(LinearEuler2D_t),intent(in) :: this - real(prec),intent(in) :: s(1:this%nvar) - real(prec),intent(in) :: nhat(1:2) - real(prec) :: exts(1:this%nvar) + subroutine AdditionalInit_LinearEuler2D_t(this) + implicit none + class(LinearEuler2D_t),intent(inout) :: this + ! Local + procedure(SELF_bcMethod),pointer :: bcfunc - exts(1) = s(1) ! density - exts(2) = (nhat(2)**2-nhat(1)**2)*s(2)-2.0_prec*nhat(1)*nhat(2)*s(3) ! u - exts(3) = (nhat(1)**2-nhat(2)**2)*s(3)-2.0_prec*nhat(1)*nhat(2)*s(2) ! v - exts(4) = s(4) ! p + bcfunc => hbc2d_NoNormalFlow_LinearEuler2D + call this%hyperbolicBCs%RegisterBoundaryCondition( & + SELF_BC_NONORMALFLOW,"no_normal_flow",bcfunc) - endfunction hbc2d_NoNormalFlow_LinearEuler2D_t + endsubroutine AdditionalInit_LinearEuler2D_t + + subroutine hbc2d_NoNormalFlow_LinearEuler2D(bc,mymodel) + !! No-normal-flow boundary condition for 2D linear Euler equations. + !! Reflects the velocity vector about the boundary normal while + !! preserving density and pressure. + class(BoundaryCondition),intent(in) :: bc + class(Model),intent(inout) :: mymodel + ! Local + integer :: n,i,iEl,j + real(prec) :: nhat(1:2),s(1:4) + + select type(m => mymodel) + class is(LinearEuler2D_t) + do n = 1,bc%nBoundaries + iEl = bc%elements(n) + j = bc%sides(n) + do i = 1,m%solution%interp%N+1 + nhat = m%geometry%nhat%boundary(i,j,iEl,1,1:2) + s = m%solution%boundary(i,j,iEl,1:4) + m%solution%extBoundary(i,j,iEl,1) = s(1) ! density + m%solution%extBoundary(i,j,iEl,2) = & + (nhat(2)**2-nhat(1)**2)*s(2)-2.0_prec*nhat(1)*nhat(2)*s(3) ! u + m%solution%extBoundary(i,j,iEl,3) = & + (nhat(1)**2-nhat(2)**2)*s(3)-2.0_prec*nhat(1)*nhat(2)*s(2) ! v + m%solution%extBoundary(i,j,iEl,4) = s(4) ! p + enddo + enddo + endselect + + endsubroutine hbc2d_NoNormalFlow_LinearEuler2D pure function flux2d_LinearEuler2D_t(this,s,dsdx) result(flux) class(LinearEuler2D_t),intent(in) :: this diff --git a/src/SELF_LinearShallowWater2D_t.f90 b/src/SELF_LinearShallowWater2D_t.f90 index 47fa9360d..d795f2714 100644 --- a/src/SELF_LinearShallowWater2D_t.f90 +++ b/src/SELF_LinearShallowWater2D_t.f90 @@ -28,6 +28,7 @@ module self_LinearShallowWater2D_t use self_model use self_dgmodel2d use self_mesh + use SELF_BoundaryConditions implicit none @@ -45,7 +46,6 @@ module self_LinearShallowWater2D_t procedure :: entropy_func => entropy_func_LinearShallowWater2D_t procedure :: flux2d => flux2d_LinearShallowWater2D_t procedure :: riemannflux2d => riemannflux2d_LinearShallowWater2D_t - procedure :: hbc2d_NoNormalFlow => hbc2d_NoNormalFlow_LinearShallowWater2D_t procedure :: sourcemethod => sourcemethod_LinearShallowWater2D_t ! Custom methods generic,public :: SetCoriolis => SetCoriolis_fplane_LinearShallowWater2D_t, & @@ -62,10 +62,16 @@ module self_LinearShallowWater2D_t subroutine AdditionalInit_LinearShallowWater2D_t(this) implicit none class(LinearShallowWater2D_t),intent(inout) :: this + ! Local + procedure(SELF_bcMethod),pointer :: bcfunc call this%fCori%Init(this%geometry%x%interp, & 1,this%mesh%nElem) + bcfunc => hbc2d_NoNormalFlow_LinearShallowWater2D + call this%hyperbolicBCs%RegisterBoundaryCondition( & + SELF_BC_NONORMALFLOW,"no_normal_flow",bcfunc) + endsubroutine AdditionalInit_LinearShallowWater2D_t subroutine AdditionalFree_LinearShallowWater2D_t(this) @@ -221,17 +227,34 @@ pure function riemannflux2d_LinearShallowWater2D_t(this,sL,sR,dsdx,nhat) result( endfunction riemannflux2d_LinearShallowWater2D_t - pure function hbc2d_NoNormalFlow_LinearShallowWater2D_t(this,s,nhat) result(exts) - class(LinearShallowWater2D_t),intent(in) :: this - real(prec),intent(in) :: s(1:this%nvar) - real(prec),intent(in) :: nhat(1:2) - real(prec) :: exts(1:this%nvar) - - exts(1) = (nhat(2)**2-nhat(1)**2)*s(1)-2.0_prec*nhat(1)*nhat(2)*s(2) ! u - exts(2) = (nhat(1)**2-nhat(2)**2)*s(2)-2.0_prec*nhat(1)*nhat(2)*s(1) ! v - exts(3) = s(3) ! eta - - endfunction hbc2d_NoNormalFlow_LinearShallowWater2D_t + subroutine hbc2d_NoNormalFlow_LinearShallowWater2D(bc,mymodel) + !! No-normal-flow boundary condition for 2D linear shallow water equations. + !! Reflects the velocity vector about the boundary normal while + !! preserving the free surface elevation. + class(BoundaryCondition),intent(in) :: bc + class(Model),intent(inout) :: mymodel + ! Local + integer :: n,i,iEl,j + real(prec) :: nhat(1:2),s(1:3) + + select type(m => mymodel) + class is(LinearShallowWater2D_t) + do n = 1,bc%nBoundaries + iEl = bc%elements(n) + j = bc%sides(n) + do i = 1,m%solution%interp%N+1 + nhat = m%geometry%nhat%boundary(i,j,iEl,1,1:2) + s = m%solution%boundary(i,j,iEl,1:3) + m%solution%extBoundary(i,j,iEl,1) = & + (nhat(2)**2-nhat(1)**2)*s(1)-2.0_prec*nhat(1)*nhat(2)*s(2) ! u + m%solution%extBoundary(i,j,iEl,2) = & + (nhat(1)**2-nhat(2)**2)*s(2)-2.0_prec*nhat(1)*nhat(2)*s(1) ! v + m%solution%extBoundary(i,j,iEl,3) = s(3) ! eta + enddo + enddo + endselect + + endsubroutine hbc2d_NoNormalFlow_LinearShallowWater2D subroutine sourcemethod_LinearShallowWater2D_t(this) implicit none diff --git a/src/SELF_Model.f90 b/src/SELF_Model.f90 index c3768df0c..8b94b1ea5 100644 --- a/src/SELF_Model.f90 +++ b/src/SELF_Model.f90 @@ -151,28 +151,6 @@ module SELF_Model procedure :: source2d => source2d_Model procedure :: source3d => source3d_Model - ! Boundary condition functions (hyperbolic) - procedure :: hbc1d_Prescribed => hbc1d_Prescribed_Model - procedure :: hbc1d_Radiation => hbc1d_Generic_Model - procedure :: hbc1d_NoNormalFlow => hbc1d_Generic_Model - procedure :: hbc2d_Prescribed => hbc2d_Prescribed_Model - procedure :: hbc2d_Radiation => hbc2d_Generic_Model - procedure :: hbc2d_NoNormalFlow => hbc2d_Generic_Model - procedure :: hbc3d_Prescribed => hbc3d_Prescribed_Model - procedure :: hbc3d_Radiation => hbc3d_Generic_Model - procedure :: hbc3d_NoNormalFlow => hbc3d_Generic_Model - - ! Boundary condition functions (parabolic) - procedure :: pbc1d_Prescribed => pbc1d_Prescribed_Model - procedure :: pbc1d_Radiation => pbc1d_Generic_Model - procedure :: pbc1d_NoNormalFlow => pbc1d_Generic_Model - procedure :: pbc2d_Prescribed => pbc2d_Prescribed_Model - procedure :: pbc2d_Radiation => pbc2d_Generic_Model - procedure :: pbc2d_NoNormalFlow => pbc2d_Generic_Model - procedure :: pbc3d_Prescribed => pbc3d_Prescribed_Model - procedure :: pbc3d_Radiation => pbc3d_Generic_Model - procedure :: pbc3d_NoNormalFlow => pbc3d_Generic_Model - procedure :: ReportEntropy => ReportEntropy_Model procedure :: ReportMetrics => ReportMetrics_Model procedure :: ReportUserMetrics => ReportUserMetrics_Model @@ -447,144 +425,6 @@ pure function source3d_Model(this,s,dsdx) result(source) endfunction source3d_Model - pure function hbc1d_Generic_Model(this,s,nhat) result(exts) - class(Model),intent(in) :: this - real(prec),intent(in) :: s(1:this%nvar) - real(prec),intent(in) :: nhat - real(prec) :: exts(1:this%nvar) - - exts = 0.0_prec - if(.false.) then ! suppress unused-dummy-argument warnings for default implementation - exts = s; exts(1) = exts(1)+nhat - endif - - endfunction hbc1d_Generic_Model - - pure function hbc1d_Prescribed_Model(this,x,t) result(exts) - class(Model),intent(in) :: this - real(prec),intent(in) :: x - real(prec),intent(in) :: t - real(prec) :: exts(1:this%nvar) - - exts = 0.0_prec - if(.false.) exts(1) = exts(1)+x+t ! suppress unused-dummy-argument warnings for default implementation - - endfunction hbc1d_Prescribed_Model - - pure function hbc2d_Generic_Model(this,s,nhat) result(exts) - class(Model),intent(in) :: this - real(prec),intent(in) :: s(1:this%nvar) - real(prec),intent(in) :: nhat(1:2) - real(prec) :: exts(1:this%nvar) - - exts = 0.0_prec - if(.false.) then ! suppress unused-dummy-argument warnings for default implementation - exts = s; exts(1) = exts(1)+nhat(1) - endif - - endfunction hbc2d_Generic_Model - - pure function hbc2d_Prescribed_Model(this,x,t) result(exts) - class(Model),intent(in) :: this - real(prec),intent(in) :: x(1:2) - real(prec),intent(in) :: t - real(prec) :: exts(1:this%nvar) - - exts = 0.0_prec - if(.false.) exts(1) = exts(1)+x(1)+t ! suppress unused-dummy-argument warnings for default implementation - - endfunction hbc2d_Prescribed_Model - - pure function hbc3d_Generic_Model(this,s,nhat) result(exts) - class(Model),intent(in) :: this - real(prec),intent(in) :: s(1:this%nvar) - real(prec),intent(in) :: nhat(1:3) - real(prec) :: exts(1:this%nvar) - - exts = 0.0_prec - if(.false.) then ! suppress unused-dummy-argument warnings for default implementation - exts = s; exts(1) = exts(1)+nhat(1) - endif - - endfunction hbc3d_Generic_Model - - pure function hbc3d_Prescribed_Model(this,x,t) result(exts) - class(Model),intent(in) :: this - real(prec),intent(in) :: x(1:3) - real(prec),intent(in) :: t - real(prec) :: exts(1:this%nvar) - - exts = 0.0_prec - if(.false.) exts(1) = exts(1)+x(1)+t ! suppress unused-dummy-argument warnings for default implementation - - endfunction hbc3d_Prescribed_Model - - pure function pbc1d_Generic_Model(this,dsdx,nhat) result(extDsdx) - class(Model),intent(in) :: this - real(prec),intent(in) :: dsdx(1:this%nvar) - real(prec),intent(in) :: nhat - real(prec) :: extDsdx(1:this%nvar) - - extDsdx = dsdx - if(.false.) extDsdx(1) = extDsdx(1)+nhat ! suppress unused-dummy-argument warning for default implementation - - endfunction pbc1d_Generic_Model - - pure function pbc1d_Prescribed_Model(this,x,t) result(extDsdx) - class(Model),intent(in) :: this - real(prec),intent(in) :: x - real(prec),intent(in) :: t - real(prec) :: extDsdx(1:this%nvar) - - extDsdx = 0.0_prec - if(.false.) extDsdx(1) = extDsdx(1)+x+t ! suppress unused-dummy-argument warnings for default implementation - - endfunction pbc1d_Prescribed_Model - - pure function pbc2d_Generic_Model(this,dsdx,nhat) result(extDsdx) - class(Model),intent(in) :: this - real(prec),intent(in) :: dsdx(1:this%nvar,1:2) - real(prec),intent(in) :: nhat(1:2) - real(prec) :: extDsdx(1:this%nvar,1:2) - - extDsdx = dsdx - if(.false.) extDsdx(1,1) = extDsdx(1,1)+nhat(1) ! suppress unused-dummy-argument warning for default implementation - - endfunction pbc2d_Generic_Model - - pure function pbc2d_Prescribed_Model(this,x,t) result(extDsdx) - class(Model),intent(in) :: this - real(prec),intent(in) :: x(1:2) - real(prec),intent(in) :: t - real(prec) :: extDsdx(1:this%nvar,1:2) - - extDsdx = 0.0_prec - if(.false.) extDsdx(1,1) = extDsdx(1,1)+x(1)+t ! suppress unused-dummy-argument warnings for default implementation - - endfunction pbc2d_Prescribed_Model - - pure function pbc3d_Generic_Model(this,dsdx,nhat) result(extDsdx) - class(Model),intent(in) :: this - real(prec),intent(in) :: dsdx(1:this%nvar,1:3) - real(prec),intent(in) :: nhat(1:3) - real(prec) :: extDsdx(1:this%nvar,1:3) - - extDsdx = dsdx - if(.false.) extDsdx(1,1) = extDsdx(1,1)+nhat(1) ! suppress unused-dummy-argument warning for default implementation - - endfunction pbc3d_Generic_Model - - pure function pbc3d_Prescribed_Model(this,x,t) result(extDsdx) - class(Model),intent(in) :: this - real(prec),intent(in) :: x(1:3) - real(prec),intent(in) :: t - real(prec) :: extDsdx(1:this%nvar,1:3) - - extDsdx = 0.0_prec - if(.false.) extDsdx(1,1) = extDsdx(1,1)+x(1)+t ! suppress unused-dummy-argument warnings for default implementation - - endfunction pbc3d_Prescribed_Model - subroutine SetTimeIntegrator_withChar(this,integrator) !! Sets the time integrator method, using a character input !! diff --git a/src/SELF_advection_diffusion_2d_t.f90 b/src/SELF_advection_diffusion_2d_t.f90 index 596dff4be..9185d6d2b 100644 --- a/src/SELF_advection_diffusion_2d_t.f90 +++ b/src/SELF_advection_diffusion_2d_t.f90 @@ -28,6 +28,7 @@ module self_advection_diffusion_2d_t use self_dgmodel2d use self_mesh + use SELF_BoundaryConditions implicit none @@ -37,6 +38,7 @@ module self_advection_diffusion_2d_t real(prec) :: v ! constant y-component of velocity contains + procedure :: AdditionalInit => AdditionalInit_advection_diffusion_2d_t procedure :: riemannflux2d => riemannflux2d_advection_diffusion_2d_t procedure :: flux2d => flux2d_advection_diffusion_2d_t procedure :: entropy_func => entropy_func_advection_diffusion_2d_t @@ -45,6 +47,84 @@ module self_advection_diffusion_2d_t contains + subroutine AdditionalInit_advection_diffusion_2d_t(this) + !! Register boundary conditions for the advection-diffusion model. + !! Hyperbolic: mirror (no-normal-flow) wall condition on the solution. + !! Parabolic: no-stress condition that zeros the normal gradient at + !! the boundary, giving zero diffusive flux through walls. This is + !! stable with the Bassi-Rebay (BR1) method. + implicit none + class(advection_diffusion_2d_t),intent(inout) :: this + ! Local + procedure(SELF_bcMethod),pointer :: bcfunc + + bcfunc => hbc2d_NoNormalFlow_advection_diffusion_2d_t + call this%hyperbolicBCs%RegisterBoundaryCondition( & + SELF_BC_NONORMALFLOW,"no_normal_flow",bcfunc) + + bcfunc => pbc2d_NoStress_advection_diffusion_2d_t + call this%parabolicBCs%RegisterBoundaryCondition( & + SELF_BC_NONORMALFLOW,"no_normal_flow",bcfunc) + + endsubroutine AdditionalInit_advection_diffusion_2d_t + + subroutine hbc2d_NoNormalFlow_advection_diffusion_2d_t(bc,mymodel) + !! Mirror boundary condition on the solution: sets the exterior + !! state equal to the interior state. With LLF this zeroes the + !! upwind dissipation at the wall. + class(BoundaryCondition),intent(in) :: bc + class(Model),intent(inout) :: mymodel + ! Local + integer :: n,i,iEl,s + + select type(m => mymodel) + class is(advection_diffusion_2d_t) + do n = 1,bc%nBoundaries + iEl = bc%elements(n) + s = bc%sides(n) + do i = 1,m%solution%interp%N+1 + m%solution%extBoundary(i,s,iEl,1:m%nvar) = & + m%solution%boundary(i,s,iEl,1:m%nvar) + enddo + enddo + endselect + + endsubroutine hbc2d_NoNormalFlow_advection_diffusion_2d_t + + subroutine pbc2d_NoStress_advection_diffusion_2d_t(bc,mymodel) + !! No-stress boundary condition for the BR1 diffusive flux. + !! Reflects the interior gradient so that the normal component + !! of the averaged gradient vanishes at the boundary: + !! sigma_ext = sigma_int - 2 (sigma_int . nhat) nhat + !! This gives zero diffusive flux through the wall and is + !! unconditionally stable for Bassi-Rebay. + class(BoundaryCondition),intent(in) :: bc + class(Model),intent(inout) :: mymodel + ! Local + integer :: n,i,iEl,s,ivar + real(prec) :: nhat(1:2),sigma_n + + select type(m => mymodel) + class is(advection_diffusion_2d_t) + do n = 1,bc%nBoundaries + iEl = bc%elements(n) + s = bc%sides(n) + do i = 1,m%solution%interp%N+1 + nhat(1:2) = m%geometry%nHat%boundary(i,s,iEl,1,1:2) + do ivar = 1,m%nvar + sigma_n = m%solutionGradient%boundary(i,s,iEl,ivar,1)*nhat(1)+ & + m%solutionGradient%boundary(i,s,iEl,ivar,2)*nhat(2) + m%solutionGradient%extBoundary(i,s,iEl,ivar,1) = & + m%solutionGradient%boundary(i,s,iEl,ivar,1)-2.0_prec*sigma_n*nhat(1) + m%solutionGradient%extBoundary(i,s,iEl,ivar,2) = & + m%solutionGradient%boundary(i,s,iEl,ivar,2)-2.0_prec*sigma_n*nhat(2) + enddo + enddo + enddo + endselect + + endsubroutine pbc2d_NoStress_advection_diffusion_2d_t + pure function entropy_func_advection_diffusion_2d_t(this,s) result(e) class(advection_diffusion_2d_t),intent(in) :: this real(prec),intent(in) :: s(1:this%solution%nvar) diff --git a/src/SELF_advection_diffusion_3d_t.f90 b/src/SELF_advection_diffusion_3d_t.f90 index 92dd69244..bbaecf676 100644 --- a/src/SELF_advection_diffusion_3d_t.f90 +++ b/src/SELF_advection_diffusion_3d_t.f90 @@ -29,6 +29,7 @@ module self_advection_diffusion_3d_t use self_model use self_dgmodel3d use self_mesh + use SELF_BoundaryConditions implicit none @@ -40,6 +41,7 @@ module self_advection_diffusion_3d_t contains + procedure :: AdditionalInit => AdditionalInit_advection_diffusion_3d_t procedure :: riemannflux3d => riemannflux3d_advection_diffusion_3d_t procedure :: flux3d => flux3d_advection_diffusion_3d_t procedure :: entropy_func => entropy_func_advection_diffusion_3d_t @@ -48,6 +50,90 @@ module self_advection_diffusion_3d_t contains + subroutine AdditionalInit_advection_diffusion_3d_t(this) + !! Register boundary conditions for the advection-diffusion model. + !! Hyperbolic: mirror (no-normal-flow) wall condition on the solution. + !! Parabolic: no-stress condition that zeros the normal gradient at + !! the boundary, giving zero diffusive flux through walls. This is + !! stable with the Bassi-Rebay (BR1) method. + implicit none + class(advection_diffusion_3d_t),intent(inout) :: this + ! Local + procedure(SELF_bcMethod),pointer :: bcfunc + + bcfunc => hbc3d_NoNormalFlow_advection_diffusion_3d_t + call this%hyperbolicBCs%RegisterBoundaryCondition( & + SELF_BC_NONORMALFLOW,"no_normal_flow",bcfunc) + + bcfunc => pbc3d_NoStress_advection_diffusion_3d_t + call this%parabolicBCs%RegisterBoundaryCondition( & + SELF_BC_NONORMALFLOW,"no_normal_flow",bcfunc) + + endsubroutine AdditionalInit_advection_diffusion_3d_t + + subroutine hbc3d_NoNormalFlow_advection_diffusion_3d_t(bc,mymodel) + !! Mirror boundary condition on the solution: sets the exterior + !! state equal to the interior state. With LLF this zeroes the + !! upwind dissipation at the wall. + class(BoundaryCondition),intent(in) :: bc + class(Model),intent(inout) :: mymodel + ! Local + integer :: n,i,j,iEl,s + + select type(m => mymodel) + class is(advection_diffusion_3d_t) + do n = 1,bc%nBoundaries + iEl = bc%elements(n) + s = bc%sides(n) + do j = 1,m%solution%interp%N+1 + do i = 1,m%solution%interp%N+1 + m%solution%extBoundary(i,j,s,iEl,1:m%nvar) = & + m%solution%boundary(i,j,s,iEl,1:m%nvar) + enddo + enddo + enddo + endselect + + endsubroutine hbc3d_NoNormalFlow_advection_diffusion_3d_t + + subroutine pbc3d_NoStress_advection_diffusion_3d_t(bc,mymodel) + !! No-stress boundary condition for the BR1 diffusive flux. + !! Reflects the interior gradient so that the normal component + !! of the averaged gradient vanishes at the boundary: + !! sigma_ext = sigma_int - 2 (sigma_int . nhat) nhat + !! This gives zero diffusive flux through the wall and is + !! unconditionally stable for Bassi-Rebay. + class(BoundaryCondition),intent(in) :: bc + class(Model),intent(inout) :: mymodel + ! Local + integer :: n,i,j,iEl,s,ivar,idir + real(prec) :: nhat(1:3),sigma_n + + select type(m => mymodel) + class is(advection_diffusion_3d_t) + do n = 1,bc%nBoundaries + iEl = bc%elements(n) + s = bc%sides(n) + do j = 1,m%solution%interp%N+1 + do i = 1,m%solution%interp%N+1 + nhat(1:3) = m%geometry%nHat%boundary(i,j,s,iEl,1,1:3) + do ivar = 1,m%nvar + sigma_n = m%solutionGradient%boundary(i,j,s,iEl,ivar,1)*nhat(1)+ & + m%solutionGradient%boundary(i,j,s,iEl,ivar,2)*nhat(2)+ & + m%solutionGradient%boundary(i,j,s,iEl,ivar,3)*nhat(3) + do idir = 1,3 + m%solutionGradient%extBoundary(i,j,s,iEl,ivar,idir) = & + m%solutionGradient%boundary(i,j,s,iEl,ivar,idir)- & + 2.0_prec*sigma_n*nhat(idir) + enddo + enddo + enddo + enddo + enddo + endselect + + endsubroutine pbc3d_NoStress_advection_diffusion_3d_t + pure function entropy_func_advection_diffusion_3d_t(this,s) result(e) class(advection_diffusion_3d_t),intent(in) :: this real(prec),intent(in) :: s(1:this%solution%nvar) diff --git a/src/gpu/SELF_DGModel1D.f90 b/src/gpu/SELF_DGModel1D.f90 index 5e027b000..393c7d6e7 100644 --- a/src/gpu/SELF_DGModel1D.f90 +++ b/src/gpu/SELF_DGModel1D.f90 @@ -37,6 +37,7 @@ module SELF_DGModel1D use SELF_DGModel1D_t use SELF_GPU use SELF_GPUInterfaces + use SELF_BoundaryConditions implicit none @@ -44,6 +45,9 @@ module SELF_DGModel1D contains + procedure :: Init => Init_DGModel1D + procedure :: Free => Free_DGModel1D + procedure :: UpdateSolution => UpdateSolution_DGModel1D procedure :: CalculateEntropy => CalculateEntropy_DGModel1D @@ -172,68 +176,90 @@ subroutine CalculateEntropy_DGModel1D(this) endsubroutine CalculateEntropy_DGModel1D - subroutine setboundarycondition_DGModel1D(this) - ! Here, we use the pre-tendency method to calculate the - ! derivative of the solution using a bassi-rebay method - ! We then do a boundary interpolation and side exchange - ! on the gradient field + subroutine Init_DGModel1D(this,mesh,geometry) + !! Initialize the 1D DG model, then upload BC element/side arrays to GPU. implicit none - class(DGModel1D),intent(inout) :: this - ! local - integer :: N,nelem - real(prec) :: x - - call gpuCheck(hipMemcpy(c_loc(this%solution%boundary), & - this%solution%boundary_gpu,sizeof(this%solution%boundary), & - hipMemcpyDeviceToHost)) - - nelem = this%geometry%nelem ! number of elements in the mesh - N = this%solution%interp%N ! polynomial degree - ! left-most boundary - if(this%mesh%bcid(1) == SELF_BC_PRESCRIBED) then - - x = this%geometry%x%boundary(1,1,1) - this%solution%extBoundary(1,1,1:this%nvar) = & - this%hbc1d_Prescribed(x,this%t) - - elseif(this%mesh%bcid(1) == SELF_BC_RADIATION) then - - this%solution%extBoundary(1,1,1:this%nvar) = & - this%hbc1d_Radiation(this%solution%boundary(1,1,1:this%nvar),-1.0_prec) - - elseif(this%mesh%bcid(1) == SELF_BC_NONORMALFLOW) then - - this%solution%extBoundary(1,1,1:this%nvar) = & - this%hbc1d_NoNormalFlow(this%solution%boundary(1,1,1:this%nvar),-1.0_prec) - - else ! Periodic - - this%solution%extBoundary(1,1,1:this%nvar) = this%solution%boundary(2,nelem,1:this%nvar) - - endif + class(DGModel1D),intent(out) :: this + type(Mesh1D),intent(in),target :: mesh + type(Geometry1D),intent(in),target :: geometry + ! Local + type(BoundaryCondition),pointer :: bc + + call Init_DGModel1D_t(this,mesh,geometry) + + ! Upload hyperbolic BC element/side arrays to device + bc => this%hyperbolicBCs%head + do while(associated(bc)) + if(bc%nBoundaries > 0) then + call gpuCheck(hipMalloc(bc%elements_gpu,sizeof(bc%elements))) + call gpuCheck(hipMemcpy(bc%elements_gpu,c_loc(bc%elements), & + sizeof(bc%elements),hipMemcpyHostToDevice)) + call gpuCheck(hipMalloc(bc%sides_gpu,sizeof(bc%sides))) + call gpuCheck(hipMemcpy(bc%sides_gpu,c_loc(bc%sides), & + sizeof(bc%sides),hipMemcpyHostToDevice)) + endif + bc => bc%next + enddo - ! right-most boundary - if(this%mesh%bcid(1) == SELF_BC_PRESCRIBED) then + ! Upload parabolic BC element/side arrays to device + bc => this%parabolicBCs%head + do while(associated(bc)) + if(bc%nBoundaries > 0) then + call gpuCheck(hipMalloc(bc%elements_gpu,sizeof(bc%elements))) + call gpuCheck(hipMemcpy(bc%elements_gpu,c_loc(bc%elements), & + sizeof(bc%elements),hipMemcpyHostToDevice)) + call gpuCheck(hipMalloc(bc%sides_gpu,sizeof(bc%sides))) + call gpuCheck(hipMemcpy(bc%sides_gpu,c_loc(bc%sides), & + sizeof(bc%sides),hipMemcpyHostToDevice)) + endif + bc => bc%next + enddo - x = this%geometry%x%boundary(2,nelem,1) - this%solution%extBoundary(2,nelem,1:this%nvar) = & - this%hbc1d_Prescribed(x,this%t) + endsubroutine Init_DGModel1D - elseif(this%mesh%bcid(1) == SELF_BC_RADIATION) then + subroutine Free_DGModel1D(this) + !! Free the 1D DG model, including GPU BC arrays. + implicit none + class(DGModel1D),intent(inout) :: this + ! Local + type(BoundaryCondition),pointer :: bc + + ! Free hyperbolic BC device arrays + bc => this%hyperbolicBCs%head + do while(associated(bc)) + if(c_associated(bc%elements_gpu)) call gpuCheck(hipFree(bc%elements_gpu)) + if(c_associated(bc%sides_gpu)) call gpuCheck(hipFree(bc%sides_gpu)) + bc%elements_gpu = c_null_ptr + bc%sides_gpu = c_null_ptr + bc => bc%next + enddo - this%solution%extBoundary(2,nelem,1:this%nvar) = & - this%hbc1d_Radiation(this%solution%boundary(2,nelem,1:this%nvar),-1.0_prec) + ! Free parabolic BC device arrays + bc => this%parabolicBCs%head + do while(associated(bc)) + if(c_associated(bc%elements_gpu)) call gpuCheck(hipFree(bc%elements_gpu)) + if(c_associated(bc%sides_gpu)) call gpuCheck(hipFree(bc%sides_gpu)) + bc%elements_gpu = c_null_ptr + bc%sides_gpu = c_null_ptr + bc => bc%next + enddo - elseif(this%mesh%bcid(1) == SELF_BC_NONORMALFLOW) then + call Free_DGModel1D_t(this) - this%solution%extBoundary(2,nelem,1:this%nvar) = & - this%hbc1d_NoNormalFlow(this%solution%boundary(2,nelem,1:this%nvar),-1.0_prec) + endsubroutine Free_DGModel1D - else ! Periodic + subroutine setboundarycondition_DGModel1D(this) + !! Apply boundary conditions for the solution on GPU. + !! Syncs boundary data from device, applies host-side BC dispatch + !! (periodic defaults + registered BCs), then syncs back to device. + implicit none + class(DGModel1D),intent(inout) :: this - this%solution%extBoundary(2,nelem,1:this%nvar) = this%solution%boundary(1,1,1:this%nvar) + call gpuCheck(hipMemcpy(c_loc(this%solution%boundary), & + this%solution%boundary_gpu,sizeof(this%solution%boundary), & + hipMemcpyDeviceToHost)) - endif + call setboundarycondition_DGModel1D_t(this) call gpuCheck(hipMemcpy(this%solution%extBoundary_gpu, & c_loc(this%solution%extBoundary), & @@ -243,68 +269,17 @@ subroutine setboundarycondition_DGModel1D(this) endsubroutine setboundarycondition_DGModel1D subroutine setgradientboundarycondition_DGModel1D(this) - ! Here, we set the boundary conditions for the - ! solution and the solution gradient at the left - ! and right most boundaries. - ! - ! Here, we use periodic boundary conditions + !! Apply gradient boundary conditions on GPU. + !! Syncs gradient boundary data from device, applies host-side BC dispatch + !! (periodic defaults + registered BCs), then syncs back to device. implicit none class(DGModel1D),intent(inout) :: this - ! local - integer :: nelem - real(prec) :: x call gpuCheck(hipMemcpy(c_loc(this%solutiongradient%boundary), & this%solutiongradient%boundary_gpu,sizeof(this%solutiongradient%boundary), & hipMemcpyDeviceToHost)) - nelem = this%geometry%nelem ! number of elements in the mesh - - ! left-most boundary - if(this%mesh%bcid(1) == SELF_BC_PRESCRIBED) then - - x = this%geometry%x%boundary(1,1,1) - this%solutionGradient%extBoundary(1,1,1:this%nvar) = & - this%pbc1d_Prescribed(x,this%t) - - elseif(this%mesh%bcid(1) == SELF_BC_RADIATION) then - - this%solutionGradient%extBoundary(1,1,1:this%nvar) = & - this%pbc1d_Radiation(this%solutionGradient%boundary(1,1,1:this%nvar),-1.0_prec) - - elseif(this%mesh%bcid(1) == SELF_BC_NONORMALFLOW) then - - this%solutionGradient%extBoundary(1,1,1:this%nvar) = & - this%pbc1d_NoNormalFlow(this%solutionGradient%boundary(1,1,1:this%nvar),-1.0_prec) - - else ! Periodic - - this%solutionGradient%extBoundary(1,1,1:this%nvar) = this%solutionGradient%boundary(2,nelem,1:this%nvar) - - endif - - ! right-most boundary - if(this%mesh%bcid(1) == SELF_BC_PRESCRIBED) then - - x = this%geometry%x%boundary(2,nelem,1) - this%solutionGradient%extBoundary(2,nelem,1:this%nvar) = & - this%pbc1d_Prescribed(x,this%t) - - elseif(this%mesh%bcid(1) == SELF_BC_RADIATION) then - - this%solutionGradient%extBoundary(2,nelem,1:this%nvar) = & - this%pbc1d_Radiation(this%solutionGradient%boundary(2,nelem,1:this%nvar),-1.0_prec) - - elseif(this%mesh%bcid(1) == SELF_BC_NONORMALFLOW) then - - this%solutionGradient%extBoundary(2,nelem,1:this%nvar) = & - this%pbc1d_NoNormalFlow(this%solutionGradient%boundary(2,nelem,1:this%nvar),-1.0_prec) - - else ! Periodic - - this%solutionGradient%extBoundary(2,nelem,1:this%nvar) = this%solutionGradient%boundary(1,1,1:this%nvar) - - endif + call setgradientboundarycondition_DGModel1D_t(this) call gpuCheck(hipMemcpy(this%solutiongradient%extBoundary_gpu, & c_loc(this%solutiongradient%extBoundary), & diff --git a/src/gpu/SELF_DGModel2D.f90 b/src/gpu/SELF_DGModel2D.f90 index aa4534912..ddc0ce9ad 100644 --- a/src/gpu/SELF_DGModel2D.f90 +++ b/src/gpu/SELF_DGModel2D.f90 @@ -29,6 +29,9 @@ module SELF_DGModel2D use SELF_DGModel2D_t use SELF_GPU use SELF_GPUInterfaces + use SELF_BoundaryConditions + use SELF_Geometry_2D + use SELF_Mesh_2D implicit none @@ -36,14 +39,15 @@ module SELF_DGModel2D contains + procedure :: Init => Init_DGModel2D + procedure :: Free => Free_DGModel2D + procedure :: UpdateSolution => UpdateSolution_DGModel2D procedure :: CalculateEntropy => CalculateEntropy_DGModel2D procedure :: BoundaryFlux => BoundaryFlux_DGModel2D procedure :: FluxMethod => fluxmethod_DGModel2D procedure :: SourceMethod => sourcemethod_DGModel2D - procedure :: SetBoundaryCondition => setboundarycondition_DGModel2D - procedure :: SetGradientBoundaryCondition => setgradientboundarycondition_DGModel2D procedure :: UpdateGRK2 => UpdateGRK2_DGModel2D procedure :: UpdateGRK3 => UpdateGRK3_DGModel2D @@ -292,136 +296,77 @@ subroutine sourcemethod_DGModel2D(this) endsubroutine sourcemethod_DGModel2D - subroutine setboundarycondition_DGModel2D(this) - !! Boundary conditions for the solution are set to - !! 0 for the external state to provide radiation type - !! boundary conditions. + subroutine Init_DGModel2D(this,mesh,geometry) + !! Initialize the 2D DG model, then upload BC element/side arrays to GPU. implicit none - class(DGModel2D),intent(inout) :: this - ! local - integer :: i,iEl,j,e2,bcid - real(prec) :: nhat(1:2),x(1:2) - - call gpuCheck(hipMemcpy(c_loc(this%solution%boundary), & - this%solution%boundary_gpu,sizeof(this%solution%boundary), & - hipMemcpyDeviceToHost)) - - call gpuCheck(hipMemcpy(c_loc(this%solution%extboundary), & - this%solution%extboundary_gpu,sizeof(this%solution%extboundary), & - hipMemcpyDeviceToHost)) - - do concurrent(j=1:4,iel=1:this%mesh%nElem) - - bcid = this%mesh%sideInfo(5,j,iEl) ! Boundary Condition ID - e2 = this%mesh%sideInfo(3,j,iEl) ! Neighboring Element ID - - if(e2 == 0) then - if(bcid == SELF_BC_PRESCRIBED) then - - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - x = this%geometry%x%boundary(i,j,iEl,1,1:2) - - this%solution%extBoundary(i,j,iEl,1:this%nvar) = & - this%hbc2d_Prescribed(x,this%t) - enddo - - elseif(bcid == SELF_BC_RADIATION) then - - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) - - this%solution%extBoundary(i,j,iEl,1:this%nvar) = & - this%hbc2d_Radiation(this%solution%boundary(i,j,iEl,1:this%nvar),nhat) - enddo - - elseif(bcid == SELF_BC_NONORMALFLOW) then - - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) - - this%solution%extBoundary(i,j,iEl,1:this%nvar) = & - this%hbc2d_NoNormalFlow(this%solution%boundary(i,j,iEl,1:this%nvar),nhat) - enddo - - endif + class(DGModel2D),intent(out) :: this + type(Mesh2D),intent(in),target :: mesh + type(SEMQuad),intent(in),target :: geometry + ! Local + type(BoundaryCondition),pointer :: bc + + call Init_DGModel2D_t(this,mesh,geometry) + + ! Upload hyperbolic BC element/side arrays to device + bc => this%hyperbolicBCs%head + do while(associated(bc)) + if(bc%nBoundaries > 0) then + call gpuCheck(hipMalloc(bc%elements_gpu,sizeof(bc%elements))) + call gpuCheck(hipMemcpy(bc%elements_gpu,c_loc(bc%elements), & + sizeof(bc%elements),hipMemcpyHostToDevice)) + call gpuCheck(hipMalloc(bc%sides_gpu,sizeof(bc%sides))) + call gpuCheck(hipMemcpy(bc%sides_gpu,c_loc(bc%sides), & + sizeof(bc%sides),hipMemcpyHostToDevice)) endif - + bc => bc%next enddo - call gpuCheck(hipMemcpy(this%solution%extBoundary_gpu, & - c_loc(this%solution%extBoundary), & - sizeof(this%solution%extBoundary), & - hipMemcpyHostToDevice)) + ! Upload parabolic BC element/side arrays to device + bc => this%parabolicBCs%head + do while(associated(bc)) + if(bc%nBoundaries > 0) then + call gpuCheck(hipMalloc(bc%elements_gpu,sizeof(bc%elements))) + call gpuCheck(hipMemcpy(bc%elements_gpu,c_loc(bc%elements), & + sizeof(bc%elements),hipMemcpyHostToDevice)) + call gpuCheck(hipMalloc(bc%sides_gpu,sizeof(bc%sides))) + call gpuCheck(hipMemcpy(bc%sides_gpu,c_loc(bc%sides), & + sizeof(bc%sides),hipMemcpyHostToDevice)) + endif + bc => bc%next + enddo - endsubroutine setboundarycondition_DGModel2D + endsubroutine Init_DGModel2D - subroutine setgradientboundarycondition_DGModel2D(this) - !! Boundary conditions for the solution are set to - !! 0 for the external state to provide radiation type - !! boundary conditions. + subroutine Free_DGModel2D(this) + !! Free the 2D DG model, including GPU BC arrays. implicit none class(DGModel2D),intent(inout) :: this - ! local - integer :: i,iEl,j,e2,bcid - real(prec) :: dsdx(1:this%nvar,1:2) - real(prec) :: nhat(1:2),x(1:2) - - call gpuCheck(hipMemcpy(c_loc(this%solutiongradient%boundary), & - this%solutiongradient%boundary_gpu,sizeof(this%solutiongradient%boundary), & - hipMemcpyDeviceToHost)) - - call gpuCheck(hipMemcpy(c_loc(this%solutiongradient%extboundary), & - this%solutiongradient%extboundary_gpu,sizeof(this%solutiongradient%extboundary), & - hipMemcpyDeviceToHost)) - - do concurrent(j=1:4,iel=1:this%mesh%nElem) - - bcid = this%mesh%sideInfo(5,j,iEl) ! Boundary Condition ID - e2 = this%mesh%sideInfo(3,j,iEl) ! Neighboring Element ID - - if(e2 == 0) then - if(bcid == SELF_BC_PRESCRIBED) then - - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - x = this%geometry%x%boundary(i,j,iEl,1,1:2) - - this%solutiongradient%extBoundary(i,j,iEl,1:this%nvar,1:2) = & - this%pbc2d_Prescribed(x,this%t) - enddo - - elseif(bcid == SELF_BC_RADIATION) then - - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) - - dsdx = this%solutiongradient%boundary(i,j,iEl,1:this%nvar,1:2) - - this%solutiongradient%extBoundary(i,j,iEl,1:this%nvar,1:2) = & - this%pbc2d_Radiation(dsdx,nhat) - enddo - - elseif(bcid == SELF_BC_NONORMALFLOW) then - - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,iEl,1,1:2) - - dsdx = this%solutiongradient%boundary(i,j,iEl,1:this%nvar,1:2) - - this%solutiongradient%extBoundary(i,j,iEl,1:this%nvar,1:2) = & - this%pbc2d_NoNormalFlow(dsdx,nhat) - enddo - - endif - endif + ! Local + type(BoundaryCondition),pointer :: bc + + ! Free hyperbolic BC device arrays + bc => this%hyperbolicBCs%head + do while(associated(bc)) + if(c_associated(bc%elements_gpu)) call gpuCheck(hipFree(bc%elements_gpu)) + if(c_associated(bc%sides_gpu)) call gpuCheck(hipFree(bc%sides_gpu)) + bc%elements_gpu = c_null_ptr + bc%sides_gpu = c_null_ptr + bc => bc%next + enddo + ! Free parabolic BC device arrays + bc => this%parabolicBCs%head + do while(associated(bc)) + if(c_associated(bc%elements_gpu)) call gpuCheck(hipFree(bc%elements_gpu)) + if(c_associated(bc%sides_gpu)) call gpuCheck(hipFree(bc%sides_gpu)) + bc%elements_gpu = c_null_ptr + bc%sides_gpu = c_null_ptr + bc => bc%next enddo - call gpuCheck(hipMemcpy(this%solutiongradient%extBoundary_gpu, & - c_loc(this%solutiongradient%extBoundary), & - sizeof(this%solutiongradient%extBoundary), & - hipMemcpyHostToDevice)) + call Free_DGModel2D_t(this) - endsubroutine setgradientboundarycondition_DGModel2D + endsubroutine Free_DGModel2D subroutine CalculateTendency_DGModel2D(this) implicit none diff --git a/src/gpu/SELF_DGModel3D.f90 b/src/gpu/SELF_DGModel3D.f90 index 908f16b09..7e79f7cea 100644 --- a/src/gpu/SELF_DGModel3D.f90 +++ b/src/gpu/SELF_DGModel3D.f90 @@ -29,6 +29,9 @@ module SELF_DGModel3D use SELF_DGModel3D_t use SELF_GPU use SELF_GPUInterfaces + use SELF_BoundaryConditions + use SELF_Geometry_3D + use SELF_Mesh_3D implicit none @@ -36,14 +39,15 @@ module SELF_DGModel3D contains + procedure :: Init => Init_DGModel3D + procedure :: Free => Free_DGModel3D + procedure :: UpdateSolution => UpdateSolution_DGModel3D procedure :: CalculateEntropy => CalculateEntropy_DGModel3D procedure :: BoundaryFlux => BoundaryFlux_DGModel3D procedure :: FluxMethod => fluxmethod_DGModel3D procedure :: SourceMethod => sourcemethod_DGModel3D - procedure :: SetBoundaryCondition => setboundarycondition_DGModel3D - procedure :: SetGradientBoundaryCondition => setgradientboundarycondition_DGModel3D procedure :: UpdateGRK2 => UpdateGRK2_DGModel3D procedure :: UpdateGRK3 => UpdateGRK3_DGModel3D @@ -293,148 +297,77 @@ subroutine sourcemethod_DGModel3D(this) endsubroutine sourcemethod_DGModel3D - subroutine setboundarycondition_DGModel3D(this) - !! Boundary conditions for the solution are set to - !! 0 for the external state to provide radiation type - !! boundary conditions. + subroutine Init_DGModel3D(this,mesh,geometry) + !! Initialize the 3D DG model, then upload BC element/side arrays to GPU. implicit none - class(DGModel3D),intent(inout) :: this - ! local - integer :: i,iEl,j,k,e2,bcid - real(prec) :: nhat(1:3),x(1:3) - - call gpuCheck(hipMemcpy(c_loc(this%solution%boundary), & - this%solution%boundary_gpu,sizeof(this%solution%boundary), & - hipMemcpyDeviceToHost)) - - call gpuCheck(hipMemcpy(c_loc(this%solution%extboundary), & - this%solution%extboundary_gpu,sizeof(this%solution%extboundary), & - hipMemcpyDeviceToHost)) - - do concurrent(k=1:6,iel=1:this%mesh%nElem) - - bcid = this%mesh%sideInfo(5,k,iEl) ! Boundary Condition ID - e2 = this%mesh%sideInfo(3,k,iEl) ! Neighboring Element ID - - if(e2 == 0) then - if(bcid == SELF_BC_PRESCRIBED) then - - do j = 1,this%solution%interp%N+1 ! Loop over quadrature points - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - x = this%geometry%x%boundary(i,j,k,iEl,1,1:3) - - this%solution%extBoundary(i,j,k,iEl,1:this%nvar) = & - this%hbc3d_Prescribed(x,this%t) - enddo - enddo - - elseif(bcid == SELF_BC_RADIATION) then - - do j = 1,this%solution%interp%N+1 ! Loop over quadrature points - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) - - this%solution%extBoundary(i,j,k,iEl,1:this%nvar) = & - this%hbc3d_Radiation(this%solution%boundary(i,j,k,iEl,1:this%nvar),nhat) - enddo - enddo - - elseif(bcid == SELF_BC_NONORMALFLOW) then - - do j = 1,this%solution%interp%N+1 ! Loop over quadrature points - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) - - this%solution%extBoundary(i,j,k,iEl,1:this%nvar) = & - this%hbc3d_NoNormalFlow(this%solution%boundary(i,j,k,iEl,1:this%nvar),nhat) - enddo - enddo - - endif + class(DGModel3D),intent(out) :: this + type(Mesh3D),intent(in),target :: mesh + type(SEMHex),intent(in),target :: geometry + ! Local + type(BoundaryCondition),pointer :: bc + + call Init_DGModel3D_t(this,mesh,geometry) + + ! Upload hyperbolic BC element/side arrays to device + bc => this%hyperbolicBCs%head + do while(associated(bc)) + if(bc%nBoundaries > 0) then + call gpuCheck(hipMalloc(bc%elements_gpu,sizeof(bc%elements))) + call gpuCheck(hipMemcpy(bc%elements_gpu,c_loc(bc%elements), & + sizeof(bc%elements),hipMemcpyHostToDevice)) + call gpuCheck(hipMalloc(bc%sides_gpu,sizeof(bc%sides))) + call gpuCheck(hipMemcpy(bc%sides_gpu,c_loc(bc%sides), & + sizeof(bc%sides),hipMemcpyHostToDevice)) endif - + bc => bc%next enddo - call gpuCheck(hipMemcpy(this%solution%extBoundary_gpu, & - c_loc(this%solution%extBoundary), & - sizeof(this%solution%extBoundary), & - hipMemcpyHostToDevice)) + ! Upload parabolic BC element/side arrays to device + bc => this%parabolicBCs%head + do while(associated(bc)) + if(bc%nBoundaries > 0) then + call gpuCheck(hipMalloc(bc%elements_gpu,sizeof(bc%elements))) + call gpuCheck(hipMemcpy(bc%elements_gpu,c_loc(bc%elements), & + sizeof(bc%elements),hipMemcpyHostToDevice)) + call gpuCheck(hipMalloc(bc%sides_gpu,sizeof(bc%sides))) + call gpuCheck(hipMemcpy(bc%sides_gpu,c_loc(bc%sides), & + sizeof(bc%sides),hipMemcpyHostToDevice)) + endif + bc => bc%next + enddo - endsubroutine setboundarycondition_DGModel3D + endsubroutine Init_DGModel3D - subroutine setgradientboundarycondition_DGModel3D(this) - !! Boundary conditions for the solution are set to - !! 0 for the external state to provide radiation type - !! boundary conditions. + subroutine Free_DGModel3D(this) + !! Free the 3D DG model, including GPU BC arrays. implicit none class(DGModel3D),intent(inout) :: this - ! local - integer :: i,iEl,j,k,e2,bcid - real(prec) :: dsdx(1:this%nvar,1:3) - real(prec) :: nhat(1:3),x(1:3) - - call gpuCheck(hipMemcpy(c_loc(this%solutiongradient%boundary), & - this%solutiongradient%boundary_gpu,sizeof(this%solutiongradient%boundary), & - hipMemcpyDeviceToHost)) - - call gpuCheck(hipMemcpy(c_loc(this%solutiongradient%extboundary), & - this%solutiongradient%extboundary_gpu,sizeof(this%solutiongradient%extboundary), & - hipMemcpyDeviceToHost)) - - do concurrent(k=1:6,iel=1:this%mesh%nElem) - - bcid = this%mesh%sideInfo(5,k,iEl) ! Boundary Condition ID - e2 = this%mesh%sideInfo(3,k,iEl) ! Neighboring Element ID - - if(e2 == 0) then - if(bcid == SELF_BC_PRESCRIBED) then - - do j = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - x = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) - - this%solutiongradient%extBoundary(i,j,k,iEl,1:this%nvar,1:3) = & - this%pbc3d_Prescribed(x,this%t) - enddo - enddo - - elseif(bcid == SELF_BC_RADIATION) then - - do j = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) - - dsdx = this%solutiongradient%boundary(i,j,k,iEl,1:this%nvar,1:3) - - this%solutiongradient%extBoundary(i,j,k,iEl,1:this%nvar,1:3) = & - this%pbc3d_Radiation(dsdx,nhat) - enddo - enddo - - elseif(bcid == SELF_BC_NONORMALFLOW) then - - do j = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - do i = 1,this%solutiongradient%interp%N+1 ! Loop over quadrature points - nhat = this%geometry%nhat%boundary(i,j,k,iEl,1,1:3) - - dsdx = this%solutiongradient%boundary(i,j,k,iEl,1:this%nvar,1:3) - - this%solutiongradient%extBoundary(i,j,k,iEl,1:this%nvar,1:3) = & - this%pbc3d_NoNormalFlow(dsdx,nhat) - enddo - enddo - - endif - endif + ! Local + type(BoundaryCondition),pointer :: bc + + ! Free hyperbolic BC device arrays + bc => this%hyperbolicBCs%head + do while(associated(bc)) + if(c_associated(bc%elements_gpu)) call gpuCheck(hipFree(bc%elements_gpu)) + if(c_associated(bc%sides_gpu)) call gpuCheck(hipFree(bc%sides_gpu)) + bc%elements_gpu = c_null_ptr + bc%sides_gpu = c_null_ptr + bc => bc%next + enddo + ! Free parabolic BC device arrays + bc => this%parabolicBCs%head + do while(associated(bc)) + if(c_associated(bc%elements_gpu)) call gpuCheck(hipFree(bc%elements_gpu)) + if(c_associated(bc%sides_gpu)) call gpuCheck(hipFree(bc%sides_gpu)) + bc%elements_gpu = c_null_ptr + bc%sides_gpu = c_null_ptr + bc => bc%next enddo - call gpuCheck(hipMemcpy(this%solutiongradient%extBoundary_gpu, & - c_loc(this%solutiongradient%extBoundary), & - sizeof(this%solutiongradient%extBoundary), & - hipMemcpyHostToDevice)) + call Free_DGModel3D_t(this) - endsubroutine setgradientboundarycondition_DGModel3D + endsubroutine Free_DGModel3D subroutine CalculateTendency_DGModel3D(this) implicit none diff --git a/src/gpu/SELF_ECAdvection2D.cpp b/src/gpu/SELF_ECAdvection2D.cpp index f821cd8e7..40a68a65b 100644 --- a/src/gpu/SELF_ECAdvection2D.cpp +++ b/src/gpu/SELF_ECAdvection2D.cpp @@ -5,35 +5,40 @@ // Mirror boundary condition: extBoundary = boundary at domain faces // ============================================================ -__global__ void setboundarycondition_ecadvection2d_gpukernel( - real *extBoundary, real *boundary, int *sideInfo, int N, int nel, int nvar) +// Mirror BC kernel for 2D EC Advection +// Sets extBoundary = boundary on pre-filtered boundary faces +__global__ void hbc2d_mirror_ecadvection2d_kernel( + real *extBoundary, real *boundary, + int *elements, int *sides, + int nBoundaries, int N, int nel, int nvar) { uint32_t idof = threadIdx.x + blockIdx.x*blockDim.x; - uint32_t ndof = (N+1)*4*nel; + uint32_t total_dofs = nBoundaries * (N+1); - if (idof < ndof) { + if (idof < total_dofs) { uint32_t i = idof % (N+1); - uint32_t s1 = (idof/(N+1)) % 4; - uint32_t e1 = idof/(N+1)/4; - uint32_t e2 = sideInfo[INDEX3(2,s1,e1,5,4)]; - if (e2 == 0) { - uint32_t ivar = blockIdx.y; - extBoundary[SCB_2D_INDEX(i,s1,e1,ivar,N,nel)] = - boundary[SCB_2D_INDEX(i,s1,e1,ivar,N,nel)]; - } + uint32_t n = idof / (N+1); + uint32_t e1 = elements[n] - 1; + uint32_t s1 = sides[n] - 1; + uint32_t ivar = blockIdx.y; + extBoundary[SCB_2D_INDEX(i,s1,e1,ivar,N,nel)] = + boundary[SCB_2D_INDEX(i,s1,e1,ivar,N,nel)]; } } extern "C" { - void setboundarycondition_ecadvection2d_gpu( - real *extBoundary, real *boundary, int *sideInfo, int N, int nel, int nvar) + void hbc2d_mirror_ecadvection2d_gpu( + real *extBoundary, real *boundary, + int *elements, int *sides, + int nBoundaries, int N, int nel, int nvar) { - int ndof = (N+1)*4*nel; + int total_dofs = nBoundaries * (N+1); int threads_per_block = 256; - int nblocks_x = ndof/threads_per_block + 1; - setboundarycondition_ecadvection2d_gpukernel<<>>(extBoundary, boundary, sideInfo, N, nel, nvar); + int nblocks_x = total_dofs/threads_per_block + 1; + hbc2d_mirror_ecadvection2d_kernel<<>>(extBoundary, boundary, + elements, sides, nBoundaries, N, nel, nvar); } } diff --git a/src/gpu/SELF_ECAdvection2D.f90 b/src/gpu/SELF_ECAdvection2D.f90 index 515ca9feb..1866cb426 100644 --- a/src/gpu/SELF_ECAdvection2D.f90 +++ b/src/gpu/SELF_ECAdvection2D.f90 @@ -27,8 +27,12 @@ module SELF_ECAdvection2D use SELF_ECAdvection2D_t + use SELF_ECDGModel2D_t use SELF_GPU use SELF_GPUInterfaces + use SELF_BoundaryConditions + use SELF_Mesh_2D + use SELF_Geometry_2D use iso_c_binding implicit none @@ -37,7 +41,9 @@ module SELF_ECAdvection2D contains - procedure :: SetBoundaryCondition => SetBoundaryCondition_ECAdvection2D + procedure :: Init => Init_ECAdvection2D + procedure :: Free => Free_ECAdvection2D + procedure :: AdditionalInit => AdditionalInit_ECAdvection2D procedure :: BoundaryFlux => BoundaryFlux_ECAdvection2D procedure :: TwoPointFluxMethod => TwoPointFluxMethod_ECAdvection2D procedure :: SourceMethod => SourceMethod_ECAdvection2D @@ -45,12 +51,13 @@ module SELF_ECAdvection2D endtype ECAdvection2D interface - subroutine setboundarycondition_ecadvection2d_gpu(extboundary,boundary,sideinfo,N,nel,nvar) & - bind(c,name="setboundarycondition_ecadvection2d_gpu") + subroutine hbc2d_mirror_ecadvection2d_gpu(extboundary,boundary, & + elements,sides,nBoundaries,N,nel,nvar) & + bind(c,name="hbc2d_mirror_ecadvection2d_gpu") use iso_c_binding - type(c_ptr),value :: extboundary,boundary,sideinfo - integer(c_int),value :: N,nel,nvar - endsubroutine setboundarycondition_ecadvection2d_gpu + type(c_ptr),value :: extboundary,boundary,elements,sides + integer(c_int),value :: nBoundaries,N,nel,nvar + endsubroutine hbc2d_mirror_ecadvection2d_gpu endinterface interface @@ -77,20 +84,87 @@ subroutine twopointfluxmethod_ecadvection2d_gpu(f,s,dsdx,u,v,N,nvar,nel) & contains - subroutine SetBoundaryCondition_ECAdvection2D(this) - !! Mirror BC on GPU: extBoundary = boundary at all domain faces. + subroutine Init_ECAdvection2D(this,mesh,geometry) + !! Initialize EC Advection 2D, then upload BC element/side arrays to GPU. + implicit none + class(ECAdvection2D),intent(out) :: this + type(Mesh2D),intent(in),target :: mesh + type(SEMQuad),intent(in),target :: geometry + ! Local + type(BoundaryCondition),pointer :: bc + + call Init_ECDGModel2D_t(this,mesh,geometry) + + ! Upload hyperbolic BC element/side arrays to device + bc => this%hyperbolicBCs%head + do while(associated(bc)) + if(bc%nBoundaries > 0) then + call gpuCheck(hipMalloc(bc%elements_gpu,sizeof(bc%elements))) + call gpuCheck(hipMemcpy(bc%elements_gpu,c_loc(bc%elements), & + sizeof(bc%elements),hipMemcpyHostToDevice)) + call gpuCheck(hipMalloc(bc%sides_gpu,sizeof(bc%sides))) + call gpuCheck(hipMemcpy(bc%sides_gpu,c_loc(bc%sides), & + sizeof(bc%sides),hipMemcpyHostToDevice)) + endif + bc => bc%next + enddo + + endsubroutine Init_ECAdvection2D + + subroutine Free_ECAdvection2D(this) + !! Free EC Advection 2D, including GPU BC arrays. implicit none class(ECAdvection2D),intent(inout) :: this + ! Local + type(BoundaryCondition),pointer :: bc - call setboundarycondition_ecadvection2d_gpu( & - this%solution%extboundary_gpu, & - this%solution%boundary_gpu, & - this%mesh%sideinfo_gpu, & - this%solution%interp%N, & - this%solution%nelem, & - this%solution%nvar) + bc => this%hyperbolicBCs%head + do while(associated(bc)) + if(c_associated(bc%elements_gpu)) call gpuCheck(hipFree(bc%elements_gpu)) + if(c_associated(bc%sides_gpu)) call gpuCheck(hipFree(bc%sides_gpu)) + bc%elements_gpu = c_null_ptr + bc%sides_gpu = c_null_ptr + bc => bc%next + enddo + + call Free_ECDGModel2D_t(this) + + endsubroutine Free_ECAdvection2D - endsubroutine SetBoundaryCondition_ECAdvection2D + subroutine AdditionalInit_ECAdvection2D(this) + implicit none + class(ECAdvection2D),intent(inout) :: this + ! Local + procedure(SELF_bcMethod),pointer :: bcfunc + + ! Call parent _t AdditionalInit (registers CPU mirror BC) + call AdditionalInit_ECAdvection2D_t(this) + + ! Re-register with GPU-accelerated version + bcfunc => hbc2d_Mirror_ECAdvection2D_GPU_wrapper + call this%hyperbolicBCs%RegisterBoundaryCondition( & + SELF_BC_NONORMALFLOW,"no_normal_flow",bcfunc) + + endsubroutine AdditionalInit_ECAdvection2D + + subroutine hbc2d_Mirror_ECAdvection2D_GPU_wrapper(bc,mymodel) + !! GPU-accelerated mirror BC for 2D EC Advection. + class(BoundaryCondition),intent(in) :: bc + class(Model),intent(inout) :: mymodel + + select type(m => mymodel) + class is(ECAdvection2D) + if(bc%nBoundaries > 0) then + call hbc2d_mirror_ecadvection2d_gpu( & + m%solution%extBoundary_gpu, & + m%solution%boundary_gpu, & + bc%elements_gpu,bc%sides_gpu, & + bc%nBoundaries,m%solution%interp%N, & + m%solution%nElem,m%solution%nvar) + endif + endselect + + endsubroutine hbc2d_Mirror_ECAdvection2D_GPU_wrapper subroutine BoundaryFlux_ECAdvection2D(this) !! LLF Riemann flux on GPU — fully device-resident. diff --git a/src/gpu/SELF_ECAdvection3D.cpp b/src/gpu/SELF_ECAdvection3D.cpp index 47a3c6174..62c7204b5 100644 --- a/src/gpu/SELF_ECAdvection3D.cpp +++ b/src/gpu/SELF_ECAdvection3D.cpp @@ -6,33 +6,43 @@ // 3D: 6 sides, boundary arrays indexed (i,j,side,iel,ivar) // ============================================================ -__global__ void setboundarycondition_ecadvection3d_gpukernel( - real *extBoundary, real *boundary, int *sideInfo, int N, int nel, int nvar) +// Mirror BC kernel for 3D EC Advection +// Sets extBoundary = boundary on pre-filtered boundary faces +__global__ void hbc3d_mirror_ecadvection3d_kernel( + real *extBoundary, real *boundary, + int *elements, int *sides, + int nBoundaries, int N, int nel, int nvar) { uint32_t idof = threadIdx.x + blockIdx.x*blockDim.x; - uint32_t ndof = (N+1)*(N+1)*6*nel; - - if (idof < ndof) { - uint32_t s1 = (idof/(N+1)/(N+1)) % 6; - uint32_t e1 = idof/(N+1)/(N+1)/6; - uint32_t e2 = sideInfo[INDEX3(2,s1,e1,5,6)]; - if (e2 == 0) { - uint32_t ivar = blockIdx.y; - extBoundary[idof + ndof*ivar] = boundary[idof + ndof*ivar]; - } + uint32_t dofs_per_face = (N+1)*(N+1); + uint32_t total_dofs = nBoundaries * dofs_per_face; + + if (idof < total_dofs) { + uint32_t i = idof % (N+1); + uint32_t j = (idof / (N+1)) % (N+1); + uint32_t n = idof / dofs_per_face; + uint32_t e1 = elements[n] - 1; + uint32_t s1 = sides[n] - 1; + uint32_t ivar = blockIdx.y; + extBoundary[SCB_3D_INDEX(i,j,s1,e1,ivar,N,nel)] = + boundary[SCB_3D_INDEX(i,j,s1,e1,ivar,N,nel)]; } } extern "C" { - void setboundarycondition_ecadvection3d_gpu( - real *extBoundary, real *boundary, int *sideInfo, int N, int nel, int nvar) + void hbc3d_mirror_ecadvection3d_gpu( + real *extBoundary, real *boundary, + int *elements, int *sides, + int nBoundaries, int N, int nel, int nvar) { - int ndof = (N+1)*(N+1)*6*nel; + int dofs_per_face = (N+1)*(N+1); + int total_dofs = nBoundaries * dofs_per_face; int threads_per_block = 256; - int nblocks_x = ndof/threads_per_block + 1; - setboundarycondition_ecadvection3d_gpukernel<<>>(extBoundary, boundary, sideInfo, N, nel, nvar); + int nblocks_x = total_dofs/threads_per_block + 1; + hbc3d_mirror_ecadvection3d_kernel<<>>(extBoundary, boundary, + elements, sides, nBoundaries, N, nel, nvar); } } diff --git a/src/gpu/SELF_ECAdvection3D.f90 b/src/gpu/SELF_ECAdvection3D.f90 index ac1e4e877..03dd62178 100644 --- a/src/gpu/SELF_ECAdvection3D.f90 +++ b/src/gpu/SELF_ECAdvection3D.f90 @@ -27,8 +27,12 @@ module SELF_ECAdvection3D use SELF_ECAdvection3D_t + use SELF_ECDGModel3D_t use SELF_GPU use SELF_GPUInterfaces + use SELF_BoundaryConditions + use SELF_Mesh_3D + use SELF_Geometry_3D use iso_c_binding implicit none @@ -37,7 +41,9 @@ module SELF_ECAdvection3D contains - procedure :: SetBoundaryCondition => SetBoundaryCondition_ECAdvection3D + procedure :: Init => Init_ECAdvection3D + procedure :: Free => Free_ECAdvection3D + procedure :: AdditionalInit => AdditionalInit_ECAdvection3D procedure :: BoundaryFlux => BoundaryFlux_ECAdvection3D procedure :: TwoPointFluxMethod => TwoPointFluxMethod_ECAdvection3D procedure :: SourceMethod => SourceMethod_ECAdvection3D @@ -45,12 +51,13 @@ module SELF_ECAdvection3D endtype ECAdvection3D interface - subroutine setboundarycondition_ecadvection3d_gpu(extboundary,boundary,sideinfo,N,nel,nvar) & - bind(c,name="setboundarycondition_ecadvection3d_gpu") + subroutine hbc3d_mirror_ecadvection3d_gpu(extboundary,boundary, & + elements,sides,nBoundaries,N,nel,nvar) & + bind(c,name="hbc3d_mirror_ecadvection3d_gpu") use iso_c_binding - type(c_ptr),value :: extboundary,boundary,sideinfo - integer(c_int),value :: N,nel,nvar - endsubroutine setboundarycondition_ecadvection3d_gpu + type(c_ptr),value :: extboundary,boundary,elements,sides + integer(c_int),value :: nBoundaries,N,nel,nvar + endsubroutine hbc3d_mirror_ecadvection3d_gpu endinterface interface @@ -77,20 +84,87 @@ subroutine twopointfluxmethod_ecadvection3d_gpu(f,s,dsdx,u,v,w,N,nvar,nel) & contains - subroutine SetBoundaryCondition_ECAdvection3D(this) - !! Mirror BC on GPU: extBoundary = boundary at all domain faces. + subroutine Init_ECAdvection3D(this,mesh,geometry) + !! Initialize EC Advection 3D, then upload BC element/side arrays to GPU. + implicit none + class(ECAdvection3D),intent(out) :: this + type(Mesh3D),intent(in),target :: mesh + type(SEMHex),intent(in),target :: geometry + ! Local + type(BoundaryCondition),pointer :: bc + + call Init_ECDGModel3D_t(this,mesh,geometry) + + ! Upload hyperbolic BC element/side arrays to device + bc => this%hyperbolicBCs%head + do while(associated(bc)) + if(bc%nBoundaries > 0) then + call gpuCheck(hipMalloc(bc%elements_gpu,sizeof(bc%elements))) + call gpuCheck(hipMemcpy(bc%elements_gpu,c_loc(bc%elements), & + sizeof(bc%elements),hipMemcpyHostToDevice)) + call gpuCheck(hipMalloc(bc%sides_gpu,sizeof(bc%sides))) + call gpuCheck(hipMemcpy(bc%sides_gpu,c_loc(bc%sides), & + sizeof(bc%sides),hipMemcpyHostToDevice)) + endif + bc => bc%next + enddo + + endsubroutine Init_ECAdvection3D + + subroutine Free_ECAdvection3D(this) + !! Free EC Advection 3D, including GPU BC arrays. implicit none class(ECAdvection3D),intent(inout) :: this + ! Local + type(BoundaryCondition),pointer :: bc - call setboundarycondition_ecadvection3d_gpu( & - this%solution%extboundary_gpu, & - this%solution%boundary_gpu, & - this%mesh%sideinfo_gpu, & - this%solution%interp%N, & - this%solution%nelem, & - this%solution%nvar) + bc => this%hyperbolicBCs%head + do while(associated(bc)) + if(c_associated(bc%elements_gpu)) call gpuCheck(hipFree(bc%elements_gpu)) + if(c_associated(bc%sides_gpu)) call gpuCheck(hipFree(bc%sides_gpu)) + bc%elements_gpu = c_null_ptr + bc%sides_gpu = c_null_ptr + bc => bc%next + enddo + + call Free_ECDGModel3D_t(this) + + endsubroutine Free_ECAdvection3D - endsubroutine SetBoundaryCondition_ECAdvection3D + subroutine AdditionalInit_ECAdvection3D(this) + implicit none + class(ECAdvection3D),intent(inout) :: this + ! Local + procedure(SELF_bcMethod),pointer :: bcfunc + + ! Call parent _t AdditionalInit (registers CPU mirror BC) + call AdditionalInit_ECAdvection3D_t(this) + + ! Re-register with GPU-accelerated version + bcfunc => hbc3d_Mirror_ECAdvection3D_GPU_wrapper + call this%hyperbolicBCs%RegisterBoundaryCondition( & + SELF_BC_NONORMALFLOW,"no_normal_flow",bcfunc) + + endsubroutine AdditionalInit_ECAdvection3D + + subroutine hbc3d_Mirror_ECAdvection3D_GPU_wrapper(bc,mymodel) + !! GPU-accelerated mirror BC for 3D EC Advection. + class(BoundaryCondition),intent(in) :: bc + class(Model),intent(inout) :: mymodel + + select type(m => mymodel) + class is(ECAdvection3D) + if(bc%nBoundaries > 0) then + call hbc3d_mirror_ecadvection3d_gpu( & + m%solution%extBoundary_gpu, & + m%solution%boundary_gpu, & + bc%elements_gpu,bc%sides_gpu, & + bc%nBoundaries,m%solution%interp%N, & + m%solution%nElem,m%solution%nvar) + endif + endselect + + endsubroutine hbc3d_Mirror_ECAdvection3D_GPU_wrapper subroutine BoundaryFlux_ECAdvection3D(this) !! LLF Riemann flux on GPU — fully device-resident. diff --git a/src/gpu/SELF_LinearEuler2D.cpp b/src/gpu/SELF_LinearEuler2D.cpp index 2880ab2e3..3250e744b 100644 --- a/src/gpu/SELF_LinearEuler2D.cpp +++ b/src/gpu/SELF_LinearEuler2D.cpp @@ -110,52 +110,89 @@ extern "C" } } -__global__ void setboundarycondition_LinearEuler2D_gpukernel(real *extBoundary, real *boundary, int *sideInfo, real *nhat, int N, int nel, int nvar){ - +// ============================================================ +// No-normal-flow BC kernel for 2D Linear Euler +// Operates on pre-filtered boundary faces via elements/sides arrays +// ============================================================ +__global__ void hbc2d_nonormalflow_lineareuler2d_kernel( + real *extBoundary, real *boundary, real *nhat, + int *elements, int *sides, + int nBoundaries, int N, int nel) +{ uint32_t idof = threadIdx.x + blockIdx.x*blockDim.x; - uint32_t ndof = (N+1)*4*nel; - - if(idof < ndof){ - uint32_t i = idof % (N+1); - uint32_t s1 = (idof/(N+1)) % 4; - uint32_t e1 = idof/(N+1)/4; - uint32_t e2 = sideInfo[INDEX3(2,s1,e1,5,4)]; - uint32_t bcid = sideInfo[INDEX3(4,s1,e1,5,4)]; - if( e2 == 0){ - if( bcid == SELF_BC_NONORMALFLOW ){ - - real u = boundary[SCB_2D_INDEX(i,s1,e1,1,N,nel)]; - real v = boundary[SCB_2D_INDEX(i,s1,e1,2,N,nel)]; - real nx = nhat[VEB_2D_INDEX(i,s1,e1,0,0,N,nel,1)]; - real ny = nhat[VEB_2D_INDEX(i,s1,e1,0,1,N,nel,1)]; - extBoundary[SCB_2D_INDEX(i,s1,e1,0,N,nel)] = boundary[SCB_2D_INDEX(i,s1,e1,0,N,nel)]; // density - extBoundary[SCB_2D_INDEX(i,s1,e1,1,N,nel)] = (ny*ny-nx*nx)*u-2.0*nx*ny*v; // u - extBoundary[SCB_2D_INDEX(i,s1,e1,2,N,nel)] = (nx*nx-ny*ny)*v-2.0*nx*ny*u; //v - extBoundary[SCB_2D_INDEX(i,s1,e1,3,N,nel)] = boundary[SCB_2D_INDEX(i,s1,e1,3,N,nel)]; // pressure - - } else if ( bcid == SELF_BC_RADIATION ){ - - extBoundary[SCB_2D_INDEX(i,s1,e1,0,N,nel)] = 0.0; - extBoundary[SCB_2D_INDEX(i,s1,e1,1,N,nel)] = 0.0; - extBoundary[SCB_2D_INDEX(i,s1,e1,2,N,nel)] = 0.0; - extBoundary[SCB_2D_INDEX(i,s1,e1,3,N,nel)] = 0.0; - - } - - } + uint32_t total_dofs = nBoundaries * (N+1); + + if(idof < total_dofs){ + uint32_t i = idof % (N+1); + uint32_t n = idof / (N+1); + uint32_t e1 = elements[n] - 1; // Fortran 1-based to C 0-based + uint32_t s1 = sides[n] - 1; + + real u = boundary[SCB_2D_INDEX(i,s1,e1,1,N,nel)]; + real v = boundary[SCB_2D_INDEX(i,s1,e1,2,N,nel)]; + real nx = nhat[VEB_2D_INDEX(i,s1,e1,0,0,N,nel,1)]; + real ny = nhat[VEB_2D_INDEX(i,s1,e1,0,1,N,nel,1)]; + + extBoundary[SCB_2D_INDEX(i,s1,e1,0,N,nel)] = boundary[SCB_2D_INDEX(i,s1,e1,0,N,nel)]; // density + extBoundary[SCB_2D_INDEX(i,s1,e1,1,N,nel)] = (ny*ny-nx*nx)*u - 2.0*nx*ny*v; // u + extBoundary[SCB_2D_INDEX(i,s1,e1,2,N,nel)] = (nx*nx-ny*ny)*v - 2.0*nx*ny*u; // v + extBoundary[SCB_2D_INDEX(i,s1,e1,3,N,nel)] = boundary[SCB_2D_INDEX(i,s1,e1,3,N,nel)]; // pressure } } -extern "C" +extern "C" { - void setboundarycondition_LinearEuler2D_gpu(real *extBoundary, real *boundary, int *sideInfo, real *nhat, int N, int nel, int nvar){ + void hbc2d_nonormalflow_lineareuler2d_gpu( + real *extBoundary, real *boundary, real *nhat, + int *elements, int *sides, + int nBoundaries, int N, int nel) + { int threads_per_block = 256; - int ndof = (N+1)*4*nel; - int nblocks_x = ndof/threads_per_block +1; + int total_dofs = nBoundaries * (N+1); + int nblocks_x = total_dofs/threads_per_block + 1; + hbc2d_nonormalflow_lineareuler2d_kernel<<>>(extBoundary, boundary, nhat, + elements, sides, nBoundaries, N, nel); + } +} - dim3 nblocks(nblocks_x,1,1); - dim3 nthreads(threads_per_block,1,1); +// ============================================================ +// Radiation BC kernel for 2D Linear Euler +// Sets extBoundary = 0 on pre-filtered boundary faces +// ============================================================ +__global__ void hbc2d_radiation_lineareuler2d_kernel( + real *extBoundary, + int *elements, int *sides, + int nBoundaries, int N, int nel) +{ + uint32_t idof = threadIdx.x + blockIdx.x*blockDim.x; + uint32_t total_dofs = nBoundaries * (N+1); + + if(idof < total_dofs){ + uint32_t i = idof % (N+1); + uint32_t n = idof / (N+1); + uint32_t e1 = elements[n] - 1; + uint32_t s1 = sides[n] - 1; + + extBoundary[SCB_2D_INDEX(i,s1,e1,0,N,nel)] = 0.0; + extBoundary[SCB_2D_INDEX(i,s1,e1,1,N,nel)] = 0.0; + extBoundary[SCB_2D_INDEX(i,s1,e1,2,N,nel)] = 0.0; + extBoundary[SCB_2D_INDEX(i,s1,e1,3,N,nel)] = 0.0; + } +} - setboundarycondition_LinearEuler2D_gpukernel<<>>(extBoundary,boundary,sideInfo,nhat,N,nel,nvar); +extern "C" +{ + void hbc2d_radiation_lineareuler2d_gpu( + real *extBoundary, + int *elements, int *sides, + int nBoundaries, int N, int nel) + { + int threads_per_block = 256; + int total_dofs = nBoundaries * (N+1); + int nblocks_x = total_dofs/threads_per_block + 1; + hbc2d_radiation_lineareuler2d_kernel<<>>(extBoundary, + elements, sides, nBoundaries, N, nel); } } diff --git a/src/gpu/SELF_LinearEuler2D.f90 b/src/gpu/SELF_LinearEuler2D.f90 index 4680a233d..c09af6464 100644 --- a/src/gpu/SELF_LinearEuler2D.f90 +++ b/src/gpu/SELF_LinearEuler2D.f90 @@ -15,7 +15,7 @@ ! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from ! this software without specific prior written permission. ! -! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT ! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT ! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT ! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY @@ -27,12 +27,14 @@ module self_LinearEuler2D use self_LinearEuler2D_t + use SELF_GPU + use SELF_BoundaryConditions implicit none type,extends(LinearEuler2D_t) :: LinearEuler2D contains - procedure :: setboundarycondition => setboundarycondition_LinearEuler2D + procedure :: AdditionalInit => AdditionalInit_LinearEuler2D procedure :: boundaryflux => boundaryflux_LinearEuler2D procedure :: fluxmethod => fluxmethod_LinearEuler2D procedure :: sourcemethod => sourcemethod_LinearEuler2D @@ -40,12 +42,23 @@ module self_LinearEuler2D endtype LinearEuler2D interface - subroutine setboundarycondition_LinearEuler2D_gpu(extboundary,boundary,sideinfo,nhat,N,nel,nvar) & - bind(c,name="setboundarycondition_LinearEuler2D_gpu") + subroutine hbc2d_nonormalflow_lineareuler2d_gpu(extboundary,boundary,nhat, & + elements,sides,nBoundaries,N,nel) & + bind(c,name="hbc2d_nonormalflow_lineareuler2d_gpu") use iso_c_binding - type(c_ptr),value :: extboundary,boundary,sideinfo,nhat - integer(c_int),value :: N,nel,nvar - endsubroutine setboundarycondition_LinearEuler2D_gpu + type(c_ptr),value :: extboundary,boundary,nhat,elements,sides + integer(c_int),value :: nBoundaries,N,nel + endsubroutine hbc2d_nonormalflow_lineareuler2d_gpu + endinterface + + interface + subroutine hbc2d_radiation_lineareuler2d_gpu(extboundary, & + elements,sides,nBoundaries,N,nel) & + bind(c,name="hbc2d_radiation_lineareuler2d_gpu") + use iso_c_binding + type(c_ptr),value :: extboundary,elements,sides + integer(c_int),value :: nBoundaries,N,nel + endsubroutine hbc2d_radiation_lineareuler2d_gpu endinterface interface @@ -72,6 +85,62 @@ subroutine boundaryflux_LinearEuler2D_gpu(fb,fextb,nhat,nscale,flux,rho0,c,N,nel contains + subroutine AdditionalInit_LinearEuler2D(this) + implicit none + class(LinearEuler2D),intent(inout) :: this + ! Local + procedure(SELF_bcMethod),pointer :: bcfunc + + ! Register GPU-accelerated BC methods, overwriting CPU versions + ! from the parent _t AdditionalInit + call AdditionalInit_LinearEuler2D_t(this) + + bcfunc => hbc2d_NoNormalFlow_LinearEuler2D_GPU_wrapper + call this%hyperbolicBCs%RegisterBoundaryCondition( & + SELF_BC_NONORMALFLOW,"no_normal_flow",bcfunc) + + bcfunc => hbc2d_Radiation_LinearEuler2D_GPU_wrapper + call this%hyperbolicBCs%RegisterBoundaryCondition( & + SELF_BC_RADIATION,"radiation",bcfunc) + + endsubroutine AdditionalInit_LinearEuler2D + + subroutine hbc2d_NoNormalFlow_LinearEuler2D_GPU_wrapper(bc,mymodel) + !! GPU-accelerated no-normal-flow BC for 2D Linear Euler. + class(BoundaryCondition),intent(in) :: bc + class(Model),intent(inout) :: mymodel + + select type(m => mymodel) + class is(LinearEuler2D) + if(bc%nBoundaries > 0) then + call hbc2d_nonormalflow_lineareuler2d_gpu( & + m%solution%extBoundary_gpu, & + m%solution%boundary_gpu, & + m%geometry%nhat%boundary_gpu, & + bc%elements_gpu,bc%sides_gpu, & + bc%nBoundaries,m%solution%interp%N,m%solution%nElem) + endif + endselect + + endsubroutine hbc2d_NoNormalFlow_LinearEuler2D_GPU_wrapper + + subroutine hbc2d_Radiation_LinearEuler2D_GPU_wrapper(bc,mymodel) + !! GPU-accelerated radiation BC for 2D Linear Euler. + class(BoundaryCondition),intent(in) :: bc + class(Model),intent(inout) :: mymodel + + select type(m => mymodel) + class is(LinearEuler2D) + if(bc%nBoundaries > 0) then + call hbc2d_radiation_lineareuler2d_gpu( & + m%solution%extBoundary_gpu, & + bc%elements_gpu,bc%sides_gpu, & + bc%nBoundaries,m%solution%interp%N,m%solution%nElem) + endif + endselect + + endsubroutine hbc2d_Radiation_LinearEuler2D_GPU_wrapper + subroutine sourcemethod_LinearEuler2D(this) implicit none class(LinearEuler2D),intent(inout) :: this @@ -108,54 +177,4 @@ subroutine fluxmethod_LinearEuler2D(this) endsubroutine fluxmethod_LinearEuler2D - subroutine setboundarycondition_LinearEuler2D(this) - !! Boundary conditions are set to periodic boundary conditions - implicit none - class(LinearEuler2D),intent(inout) :: this - ! local - integer :: i,iEl,j,e2,bcid - real(prec) :: x(1:2) - - if(this%prescribed_bcs_enabled) then - call gpuCheck(hipMemcpy(c_loc(this%solution%extboundary), & - this%solution%extboundary_gpu,sizeof(this%solution%extboundary), & - hipMemcpyDeviceToHost)) - - ! Prescribed boundaries are still done on the GPU - do iEl = 1,this%solution%nElem ! Loop over all elements - do j = 1,4 ! Loop over all sides - - bcid = this%mesh%sideInfo(5,j,iEl) ! Boundary Condition ID - e2 = this%mesh%sideInfo(3,j,iEl) ! Neighboring Element ID - - if(e2 == 0) then - if(bcid == SELF_BC_PRESCRIBED) then - - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - x = this%geometry%x%boundary(i,j,iEl,1,1:2) - - this%solution%extBoundary(i,j,iEl,1:this%nvar) = & - this%hbc2d_Prescribed(x,this%t) - enddo - - endif - endif - - enddo - enddo - - call gpuCheck(hipMemcpy(this%solution%extBoundary_gpu, & - c_loc(this%solution%extBoundary), & - sizeof(this%solution%extBoundary), & - hipMemcpyHostToDevice)) - endif - call setboundarycondition_LinearEuler2D_gpu(this%solution%extboundary_gpu, & - this%solution%boundary_gpu, & - this%mesh%sideInfo_gpu, & - this%geometry%nhat%boundary_gpu, & - this%solution%interp%N, & - this%solution%nelem,this%solution%nvar) - - endsubroutine setboundarycondition_LinearEuler2D - endmodule self_LinearEuler2D diff --git a/src/gpu/SELF_LinearEuler3D.cpp b/src/gpu/SELF_LinearEuler3D.cpp index c71c4b58f..0d61687e0 100644 --- a/src/gpu/SELF_LinearEuler3D.cpp +++ b/src/gpu/SELF_LinearEuler3D.cpp @@ -122,55 +122,45 @@ extern "C" } } -__global__ void setboundarycondition_LinearEuler3D_gpukernel(real *extBoundary, real *boundary, int *sideInfo, real *nhat, int N, int nel){ - +// Radiation BC kernel for 3D Linear Euler +// Sets extBoundary = 0 on pre-filtered boundary faces +__global__ void hbc3d_radiation_lineareuler3d_kernel( + real *extBoundary, + int *elements, int *sides, + int nBoundaries, int N, int nel) +{ uint32_t idof = threadIdx.x + blockIdx.x*blockDim.x; - uint32_t ndof = (N+1)*(N+1)*6*nel; - - if(idof < ndof){ - uint32_t i = idof % (N+1); - uint32_t j = (idof/(N+1)) % (N+1); - uint32_t s1 = (idof/(N+1)/(N+1)) % 6; - uint32_t e1 = idof/(N+1)/(N+1)/6; - uint32_t e2 = sideInfo[INDEX3(2,s1,e1,5,4)]; - uint32_t bcid = sideInfo[INDEX3(4,s1,e1,5,4)]; - if( e2 == 0){ - // if( bcid == SELF_BC_NONORMALFLOW ){ - - // real u = boundary[SCB_3D_INDEX(i,s1,e1,1,N,nel)]; - // real v = boundary[SCB_3D_INDEX(i,s1,e1,2,N,nel)]; - // real nx = nhat[VEB_3D_INDEX(i,s1,e1,0,0,N,nel,1)]; - // real ny = nhat[VEB_3D_INDEX(i,s1,e1,0,1,N,nel,1)]; - // extBoundary[SCB_3D_INDEX(i,s1,e1,0,N,nel)] = boundary[SCB_3D_INDEX(i,s1,e1,0,N,nel)]; // density - // extBoundary[SCB_3D_INDEX(i,s1,e1,1,N,nel)] = (ny*ny-nx*nx)*u-2.0*nx*ny*v; // u - // extBoundary[SCB_3D_INDEX(i,s1,e1,2,N,nel)] = (nx*nx-ny*ny)*v-2.0*nx*ny*u; //v - // extBoundary[SCB_3D_INDEX(i,s1,e1,3,N,nel)] = boundary[SCB_3D_INDEX(i,s1,e1,3,N,nel)]; // pressure - - // } else - if ( bcid == SELF_BC_RADIATION ){ - - extBoundary[SCB_3D_INDEX(i,j,s1,e1,0,N,nel)] = 0.0; - extBoundary[SCB_3D_INDEX(i,j,s1,e1,1,N,nel)] = 0.0; - extBoundary[SCB_3D_INDEX(i,j,s1,e1,2,N,nel)] = 0.0; - extBoundary[SCB_3D_INDEX(i,j,s1,e1,3,N,nel)] = 0.0; - extBoundary[SCB_3D_INDEX(i,j,s1,e1,4,N,nel)] = 0.0; - - } - - } + uint32_t dofs_per_face = (N+1)*(N+1); + uint32_t total_dofs = nBoundaries * dofs_per_face; + + if(idof < total_dofs){ + uint32_t i = idof % (N+1); + uint32_t j = (idof / (N+1)) % (N+1); + uint32_t n = idof / dofs_per_face; + uint32_t e1 = elements[n] - 1; + uint32_t s1 = sides[n] - 1; + + extBoundary[SCB_3D_INDEX(i,j,s1,e1,0,N,nel)] = 0.0; + extBoundary[SCB_3D_INDEX(i,j,s1,e1,1,N,nel)] = 0.0; + extBoundary[SCB_3D_INDEX(i,j,s1,e1,2,N,nel)] = 0.0; + extBoundary[SCB_3D_INDEX(i,j,s1,e1,3,N,nel)] = 0.0; + extBoundary[SCB_3D_INDEX(i,j,s1,e1,4,N,nel)] = 0.0; } } -extern "C" +extern "C" { - void setboundarycondition_LinearEuler3D_gpu(real *extBoundary, real *boundary, int *sideInfo, real *nhat, int N, int nel){ + void hbc3d_radiation_lineareuler3d_gpu( + real *extBoundary, + int *elements, int *sides, + int nBoundaries, int N, int nel) + { int threads_per_block = 256; - int ndof = (N+1)*(N+1)*6*nel; - int nblocks_x = ndof/threads_per_block +1; - - dim3 nblocks(nblocks_x,1,1); - dim3 nthreads(threads_per_block,1,1); - - setboundarycondition_LinearEuler3D_gpukernel<<>>(extBoundary,boundary,sideInfo,nhat,N,nel); + int dofs_per_face = (N+1)*(N+1); + int total_dofs = nBoundaries * dofs_per_face; + int nblocks_x = total_dofs/threads_per_block + 1; + hbc3d_radiation_lineareuler3d_kernel<<>>(extBoundary, + elements, sides, nBoundaries, N, nel); } } diff --git a/src/gpu/SELF_LinearEuler3D.f90 b/src/gpu/SELF_LinearEuler3D.f90 index 76190a383..4da12fb41 100644 --- a/src/gpu/SELF_LinearEuler3D.f90 +++ b/src/gpu/SELF_LinearEuler3D.f90 @@ -15,7 +15,7 @@ ! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from ! this software without specific prior written permission. ! -! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT ! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT ! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT ! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY @@ -27,24 +27,27 @@ module self_LinearEuler3D use self_LinearEuler3D_t + use SELF_GPU + use SELF_BoundaryConditions implicit none type,extends(LinearEuler3D_t) :: LinearEuler3D contains - procedure :: setboundarycondition => setboundarycondition_LinearEuler3D + procedure :: AdditionalInit => AdditionalInit_LinearEuler3D procedure :: boundaryflux => boundaryflux_LinearEuler3D procedure :: fluxmethod => fluxmethod_LinearEuler3D endtype LinearEuler3D interface - subroutine setboundarycondition_LinearEuler3D_gpu(extboundary,boundary,sideinfo,nhat,N,nel) & - bind(c,name="setboundarycondition_LinearEuler3D_gpu") + subroutine hbc3d_radiation_lineareuler3d_gpu(extboundary, & + elements,sides,nBoundaries,N,nel) & + bind(c,name="hbc3d_radiation_lineareuler3d_gpu") use iso_c_binding - type(c_ptr),value :: extboundary,boundary,sideinfo,nhat - integer(c_int),value :: N,nel - endsubroutine setboundarycondition_LinearEuler3D_gpu + type(c_ptr),value :: extboundary,elements,sides + integer(c_int),value :: nBoundaries,N,nel + endsubroutine hbc3d_radiation_lineareuler3d_gpu endinterface interface @@ -71,6 +74,34 @@ subroutine boundaryflux_LinearEuler3D_gpu(fb,fextb,nhat,nscale,flux,rho0,c,N,nel contains + subroutine AdditionalInit_LinearEuler3D(this) + implicit none + class(LinearEuler3D),intent(inout) :: this + ! Local + procedure(SELF_bcMethod),pointer :: bcfunc + + bcfunc => hbc3d_Radiation_LinearEuler3D_GPU_wrapper + call this%hyperbolicBCs%RegisterBoundaryCondition( & + SELF_BC_RADIATION,"radiation",bcfunc) + + endsubroutine AdditionalInit_LinearEuler3D + + subroutine hbc3d_Radiation_LinearEuler3D_GPU_wrapper(bc,mymodel) + class(BoundaryCondition),intent(in) :: bc + class(Model),intent(inout) :: mymodel + + select type(m => mymodel) + class is(LinearEuler3D) + if(bc%nBoundaries > 0) then + call hbc3d_radiation_lineareuler3d_gpu( & + m%solution%extBoundary_gpu, & + bc%elements_gpu,bc%sides_gpu, & + bc%nBoundaries,m%solution%interp%N,m%solution%nElem) + endif + endselect + + endsubroutine hbc3d_Radiation_LinearEuler3D_GPU_wrapper + subroutine boundaryflux_LinearEuler3D(this) implicit none class(LinearEuler3D),intent(inout) :: this @@ -96,56 +127,4 @@ subroutine fluxmethod_LinearEuler3D(this) endsubroutine fluxmethod_LinearEuler3D - subroutine setboundarycondition_LinearEuler3D(this) - !! Boundary conditions are set to periodic boundary conditions - implicit none - class(LinearEuler3D),intent(inout) :: this - ! local - integer :: i,iEl,j,k,e2,bcid - real(prec) :: x(1:3) - - if(this%prescribed_bcs_enabled) then - call gpuCheck(hipMemcpy(c_loc(this%solution%extboundary), & - this%solution%extboundary_gpu,sizeof(this%solution%extboundary), & - hipMemcpyDeviceToHost)) - - ! Prescribed boundaries are still done on the CPU - do iEl = 1,this%solution%nElem ! Loop over all elements - do k = 1,6 ! Loop over all sides - - bcid = this%mesh%sideInfo(5,j,iEl) ! Boundary Condition ID - e2 = this%mesh%sideInfo(3,j,iEl) ! Neighboring Element ID - - if(e2 == 0) then - if(bcid == SELF_BC_PRESCRIBED) then - - do j = 1,this%solution%interp%N+1 ! Loop over quadrature points - do i = 1,this%solution%interp%N+1 ! Loop over quadrature points - x = this%geometry%x%boundary(i,j,k,iEl,1,1:3) - - this%solution%extBoundary(i,j,k,iEl,1:this%nvar) = & - this%hbc3D_Prescribed(x,this%t) - enddo - enddo - - endif - endif - - enddo - enddo - - call gpuCheck(hipMemcpy(this%solution%extBoundary_gpu, & - c_loc(this%solution%extBoundary), & - sizeof(this%solution%extBoundary), & - hipMemcpyHostToDevice)) - endif - call setboundarycondition_LinearEuler3D_gpu(this%solution%extboundary_gpu, & - this%solution%boundary_gpu, & - this%mesh%sideInfo_gpu, & - this%geometry%nhat%boundary_gpu, & - this%solution%interp%N, & - this%solution%nelem) - - endsubroutine setboundarycondition_LinearEuler3D - endmodule self_LinearEuler3D diff --git a/src/gpu/SELF_LinearShallowWater2D.cpp b/src/gpu/SELF_LinearShallowWater2D.cpp index 7d68dbeab..b44baab63 100644 --- a/src/gpu/SELF_LinearShallowWater2D.cpp +++ b/src/gpu/SELF_LinearShallowWater2D.cpp @@ -100,47 +100,83 @@ extern "C" } } -__global__ void setboundarycondition_LinearShallowWater2D_gpukernel(real *extBoundary, real *boundary, int *sideInfo, real *nhat, int N, int nEl, int nvar){ - uint32_t idof = threadIdx.x + blockIdx.x*blockDim.x; - uint32_t ndof = (N+1)*4*nEl; - - if(idof < ndof){ - uint32_t i = idof % (N+1); - uint32_t s1 = (idof/(N+1)) % 4; - uint32_t e1 = idof/(N+1)/4; - uint32_t e2 = sideInfo[INDEX3(2,s1,e1,5,4)]; - uint32_t bcid = sideInfo[INDEX3(4,s1,e1,5,4)]; - if( e2 == 0){ - if( bcid == SELF_BC_NONORMALFLOW){ - real u = boundary[SCB_2D_INDEX(i,s1,e1,0,N,nEl)]; - real v = boundary[SCB_2D_INDEX(i,s1,e1,1,N,nEl)]; - real eta = boundary[SCB_2D_INDEX(i,s1,e1,2,N,nEl)]; - real nx = nhat[VEB_2D_INDEX(i,s1,e1,0,0,N,nEl,1)]; - real ny = nhat[VEB_2D_INDEX(i,s1,e1,0,1,N,nEl,1)]; - - extBoundary[SCB_2D_INDEX(i,s1,e1,0,N,nEl)] = (ny * ny - nx * nx) * u - 2 * nx * ny * v; - extBoundary[SCB_2D_INDEX(i,s1,e1,1,N,nEl)] = (nx * nx - ny * ny) * v - 2 * nx * ny * u; - extBoundary[SCB_2D_INDEX(i,s1,e1,2,N,nEl)] = eta; - } else if ( bcid == SELF_BC_RADIATION){ - extBoundary[SCB_2D_INDEX(i,s1,e1,0,N,nEl)] = 0.0; - extBoundary[SCB_2D_INDEX(i,s1,e1,1,N,nEl)] = 0.0; - extBoundary[SCB_2D_INDEX(i,s1,e1,2,N,nEl)] = 0.0; - } - } - } +// No-normal-flow BC kernel for 2D Linear Shallow Water +__global__ void hbc2d_nonormalflow_linearshallowwater2d_kernel( + real *extBoundary, real *boundary, real *nhat, + int *elements, int *sides, + int nBoundaries, int N, int nel) +{ + uint32_t idof = threadIdx.x + blockIdx.x*blockDim.x; + uint32_t total_dofs = nBoundaries * (N+1); + + if(idof < total_dofs){ + uint32_t i = idof % (N+1); + uint32_t n = idof / (N+1); + uint32_t e1 = elements[n] - 1; + uint32_t s1 = sides[n] - 1; + + real u = boundary[SCB_2D_INDEX(i,s1,e1,0,N,nel)]; + real v = boundary[SCB_2D_INDEX(i,s1,e1,1,N,nel)]; + real eta = boundary[SCB_2D_INDEX(i,s1,e1,2,N,nel)]; + real nx = nhat[VEB_2D_INDEX(i,s1,e1,0,0,N,nel,1)]; + real ny = nhat[VEB_2D_INDEX(i,s1,e1,0,1,N,nel,1)]; + + extBoundary[SCB_2D_INDEX(i,s1,e1,0,N,nel)] = (ny*ny - nx*nx)*u - 2.0*nx*ny*v; + extBoundary[SCB_2D_INDEX(i,s1,e1,1,N,nel)] = (nx*nx - ny*ny)*v - 2.0*nx*ny*u; + extBoundary[SCB_2D_INDEX(i,s1,e1,2,N,nel)] = eta; + } } -extern "C" +extern "C" { - void setboundarycondition_LinearShallowWater2D_gpu(real *extBoundary, real *boundary, int *sideInfo, real *nhat, int N, int nel, int nvar){ + void hbc2d_nonormalflow_linearshallowwater2d_gpu( + real *extBoundary, real *boundary, real *nhat, + int *elements, int *sides, + int nBoundaries, int N, int nel) + { int threads_per_block = 256; - int ndof = (N+1)*4*nel; - int nblocks_x = ndof/threads_per_block +1; + int total_dofs = nBoundaries * (N+1); + int nblocks_x = total_dofs/threads_per_block + 1; + hbc2d_nonormalflow_linearshallowwater2d_kernel<<>>(extBoundary, boundary, nhat, + elements, sides, nBoundaries, N, nel); + } +} - dim3 nblocks(nblocks_x,1,1); - dim3 nthreads(threads_per_block,1,1); +// Radiation BC kernel for 2D Linear Shallow Water +__global__ void hbc2d_radiation_linearshallowwater2d_kernel( + real *extBoundary, + int *elements, int *sides, + int nBoundaries, int N, int nel) +{ + uint32_t idof = threadIdx.x + blockIdx.x*blockDim.x; + uint32_t total_dofs = nBoundaries * (N+1); + + if(idof < total_dofs){ + uint32_t i = idof % (N+1); + uint32_t n = idof / (N+1); + uint32_t e1 = elements[n] - 1; + uint32_t s1 = sides[n] - 1; - setboundarycondition_LinearShallowWater2D_gpukernel<<>>(extBoundary,boundary,sideInfo,nhat,N,nel,nvar); + extBoundary[SCB_2D_INDEX(i,s1,e1,0,N,nel)] = 0.0; + extBoundary[SCB_2D_INDEX(i,s1,e1,1,N,nel)] = 0.0; + extBoundary[SCB_2D_INDEX(i,s1,e1,2,N,nel)] = 0.0; + } +} + +extern "C" +{ + void hbc2d_radiation_linearshallowwater2d_gpu( + real *extBoundary, + int *elements, int *sides, + int nBoundaries, int N, int nel) + { + int threads_per_block = 256; + int total_dofs = nBoundaries * (N+1); + int nblocks_x = total_dofs/threads_per_block + 1; + hbc2d_radiation_linearshallowwater2d_kernel<<>>(extBoundary, + elements, sides, nBoundaries, N, nel); } } diff --git a/src/gpu/SELF_LinearShallowWater2D.f90 b/src/gpu/SELF_LinearShallowWater2D.f90 index 9ae072d4f..3eef6e135 100644 --- a/src/gpu/SELF_LinearShallowWater2D.f90 +++ b/src/gpu/SELF_LinearShallowWater2D.f90 @@ -15,7 +15,7 @@ ! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from ! this software without specific prior written permission. ! -! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT ! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT ! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT ! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY @@ -27,12 +27,14 @@ module self_LinearShallowWater2D use self_LinearShallowWater2D_t + use SELF_GPU + use SELF_BoundaryConditions implicit none type,extends(LinearShallowWater2D_t) :: LinearShallowWater2D contains - procedure :: setboundarycondition => setboundarycondition_LinearShallowWater2D + procedure :: AdditionalInit => AdditionalInit_LinearShallowWater2D procedure :: boundaryflux => boundaryflux_LinearShallowWater2D procedure :: fluxmethod => fluxmethod_LinearShallowWater2D procedure :: sourcemethod => sourcemethod_LinearShallowWater2D @@ -40,12 +42,23 @@ module self_LinearShallowWater2D endtype LinearShallowWater2D interface - subroutine setboundarycondition_LinearShallowWater2D_gpu(extboundary,boundary,sideinfo,nhat,N,nel,nvar) & - bind(c,name="setboundarycondition_LinearShallowWater2D_gpu") + subroutine hbc2d_nonormalflow_linearshallowwater2d_gpu(extboundary,boundary,nhat, & + elements,sides,nBoundaries,N,nel) & + bind(c,name="hbc2d_nonormalflow_linearshallowwater2d_gpu") use iso_c_binding - type(c_ptr),value :: extboundary,boundary,sideinfo,nhat - integer(c_int),value :: N,nel,nvar - endsubroutine setboundarycondition_LinearShallowWater2D_gpu + type(c_ptr),value :: extboundary,boundary,nhat,elements,sides + integer(c_int),value :: nBoundaries,N,nel + endsubroutine hbc2d_nonormalflow_linearshallowwater2d_gpu + endinterface + + interface + subroutine hbc2d_radiation_linearshallowwater2d_gpu(extboundary, & + elements,sides,nBoundaries,N,nel) & + bind(c,name="hbc2d_radiation_linearshallowwater2d_gpu") + use iso_c_binding + type(c_ptr),value :: extboundary,elements,sides + integer(c_int),value :: nBoundaries,N,nel + endsubroutine hbc2d_radiation_linearshallowwater2d_gpu endinterface interface @@ -83,6 +96,60 @@ subroutine sourcemethod_LinearShallowWater2D_gpu(solution,source,fCori,Cd,N,nel, contains + subroutine AdditionalInit_LinearShallowWater2D(this) + implicit none + class(LinearShallowWater2D),intent(inout) :: this + ! Local + procedure(SELF_bcMethod),pointer :: bcfunc + + ! Call parent _t AdditionalInit (registers CPU BCs + initializes fCori) + call AdditionalInit_LinearShallowWater2D_t(this) + + ! Re-register with GPU-accelerated versions + bcfunc => hbc2d_NoNormalFlow_LinearShallowWater2D_GPU_wrapper + call this%hyperbolicBCs%RegisterBoundaryCondition( & + SELF_BC_NONORMALFLOW,"no_normal_flow",bcfunc) + + bcfunc => hbc2d_Radiation_LinearShallowWater2D_GPU_wrapper + call this%hyperbolicBCs%RegisterBoundaryCondition( & + SELF_BC_RADIATION,"radiation",bcfunc) + + endsubroutine AdditionalInit_LinearShallowWater2D + + subroutine hbc2d_NoNormalFlow_LinearShallowWater2D_GPU_wrapper(bc,mymodel) + class(BoundaryCondition),intent(in) :: bc + class(Model),intent(inout) :: mymodel + + select type(m => mymodel) + class is(LinearShallowWater2D) + if(bc%nBoundaries > 0) then + call hbc2d_nonormalflow_linearshallowwater2d_gpu( & + m%solution%extBoundary_gpu, & + m%solution%boundary_gpu, & + m%geometry%nhat%boundary_gpu, & + bc%elements_gpu,bc%sides_gpu, & + bc%nBoundaries,m%solution%interp%N,m%solution%nElem) + endif + endselect + + endsubroutine hbc2d_NoNormalFlow_LinearShallowWater2D_GPU_wrapper + + subroutine hbc2d_Radiation_LinearShallowWater2D_GPU_wrapper(bc,mymodel) + class(BoundaryCondition),intent(in) :: bc + class(Model),intent(inout) :: mymodel + + select type(m => mymodel) + class is(LinearShallowWater2D) + if(bc%nBoundaries > 0) then + call hbc2d_radiation_linearshallowwater2d_gpu( & + m%solution%extBoundary_gpu, & + bc%elements_gpu,bc%sides_gpu, & + bc%nBoundaries,m%solution%interp%N,m%solution%nElem) + endif + endselect + + endsubroutine hbc2d_Radiation_LinearShallowWater2D_GPU_wrapper + subroutine boundaryflux_LinearShallowWater2D(this) implicit none class(LinearShallowWater2D),intent(inout) :: this @@ -92,67 +159,20 @@ subroutine boundaryflux_LinearShallowWater2D(this) this%geometry%nhat%boundary_gpu, & this%geometry%nscale%boundary_gpu, & this%flux%boundaryNormal_gpu, & - this%g, & - this%H, & + this%g,this%H, & this%solution%interp%N, & this%solution%nelem, & this%solution%nvar) endsubroutine boundaryflux_LinearShallowWater2D - subroutine setboundarycondition_LinearShallowWater2D(this) - implicit none - class(LinearShallowWater2D),intent(inout) :: this - integer :: i,iel,j,e2,bcid - real(prec) :: x(1:2) - - if(this%prescribed_bcs_enabled) then - call gpuCheck(hipMemcpy(c_loc(this%solution%extBoundary), & - this%solution%extBoundary_gpu, & - sizeof(this%solution%extBoundary), & - hipMemcpyDeviceToHost)) - do iel = 1,this%solution%nelem - do j = 1,4 - bcid = this%mesh%sideinfo(5,j,iel) - e2 = this%mesh%sideinfo(3,j,iel) - - if(e2 == 0) then - if(bcid == SELF_BC_PRESCRIBED) then - do i = 1,this%solution%interp%N+1 - x = this%geometry%x%boundary(i,j,iel,1,1:2) - this%solution%extBoundary(i,j,iel,1:this%nvar) = & - this%hbc2d_Prescribed(x,this%t) - enddo - endif - endif - enddo - enddo - - call gpucheck(hipMemcpy(this%solution%extBoundary_gpu, & - c_loc(this%solution%extBoundary), & - sizeof(this%solution%extBoundary), & - hipMemcpyHostToDevice)) - - endif - - call setboundarycondition_LinearShallowWater2D_gpu(this%solution%extboundary_gpu, & - this%solution%boundary_gpu, & - this%mesh%sideInfo_gpu, & - this%geometry%nhat%boundary_gpu, & - this%solution%interp%N, & - this%solution%nelem, & - this%solution%nvar) - - endsubroutine setboundarycondition_LinearShallowWater2D - subroutine fluxmethod_LinearShallowWater2D(this) implicit none class(LinearShallowWater2D),intent(inout) :: this call fluxmethod_LinearShallowWater2D_gpu(this%solution%interior_gpu, & this%flux%interior_gpu, & - this%g, & - this%H, & + this%g,this%H, & this%solution%interp%N, & this%solution%nelem, & this%solution%nvar) diff --git a/src/gpu/SELF_advection_diffusion_1d.f90 b/src/gpu/SELF_advection_diffusion_1d.f90 index 126408f1a..fe1ffa4f0 100644 --- a/src/gpu/SELF_advection_diffusion_1d.f90 +++ b/src/gpu/SELF_advection_diffusion_1d.f90 @@ -35,22 +35,11 @@ module self_advection_diffusion_1d type,extends(advection_diffusion_1d_t) :: advection_diffusion_1d contains - procedure :: setboundarycondition => setboundarycondition_advection_diffusion_1d - procedure :: setgradientboundarycondition => setgradientboundarycondition_advection_diffusion_1d procedure :: boundaryflux => boundaryflux_advection_diffusion_1d procedure :: fluxmethod => fluxmethod_advection_diffusion_1d endtype advection_diffusion_1d - interface - subroutine setboundarycondition_advection_diffusion_1d_gpu(extboundary,boundary,nel,nvar) & - bind(c,name="setboundarycondition_advection_diffusion_1d_gpu") - use iso_c_binding - type(c_ptr),value :: extboundary,boundary - integer(c_int),value :: nel,nvar - endsubroutine setboundarycondition_advection_diffusion_1d_gpu - endinterface - interface subroutine fluxmethod_advection_diffusion_1d_gpu(solution,solutiongradient,flux,u,nu,ndof) & bind(c,name="fluxmethod_advection_diffusion_1d_gpu") @@ -75,26 +64,6 @@ subroutine boundaryflux_advection_diffusion_1d_gpu(fb,fextb,dfavg,flux,u,nu,ndof contains - subroutine setboundarycondition_advection_diffusion_1d(this) - !! Boundary conditions are set to periodic boundary conditions - implicit none - class(advection_diffusion_1d),intent(inout) :: this - - call setboundarycondition_advection_diffusion_1d_gpu(this%solution%extboundary_gpu, & - this%solution%boundary_gpu,this%solution%nelem,this%solution%nvar) - - endsubroutine setboundarycondition_advection_diffusion_1d - - subroutine setgradientboundarycondition_advection_diffusion_1d(this) - !! Gradient boundary conditions are set to periodic boundary conditions - implicit none - class(advection_diffusion_1d),intent(inout) :: this - - call setboundarycondition_advection_diffusion_1d_gpu(this%solutiongradient%extboundary_gpu, & - this%solutiongradient%boundary_gpu,this%solution%nelem,this%solution%nvar) - - endsubroutine setgradientboundarycondition_advection_diffusion_1d - subroutine fluxmethod_advection_diffusion_1d(this) implicit none class(advection_diffusion_1d),intent(inout) :: this diff --git a/src/gpu/SELF_advection_diffusion_2d.f90 b/src/gpu/SELF_advection_diffusion_2d.f90 index 681ddafe4..3fe82f53b 100644 --- a/src/gpu/SELF_advection_diffusion_2d.f90 +++ b/src/gpu/SELF_advection_diffusion_2d.f90 @@ -33,31 +33,11 @@ module self_advection_diffusion_2d type,extends(advection_diffusion_2d_t) :: advection_diffusion_2d contains - procedure :: setboundarycondition => setboundarycondition_advection_diffusion_2d - procedure :: setgradientboundarycondition => setgradientboundarycondition_advection_diffusion_2d procedure :: boundaryflux => boundaryflux_advection_diffusion_2d procedure :: fluxmethod => fluxmethod_advection_diffusion_2d endtype advection_diffusion_2d - interface - subroutine setboundarycondition_advection_diffusion_2d_gpu(extboundary,boundary,sideinfo,N,nel,nvar) & - bind(c,name="setboundarycondition_advection_diffusion_2d_gpu") - use iso_c_binding - type(c_ptr),value :: extboundary,boundary,sideinfo - integer(c_int),value :: N,nel,nvar - endsubroutine setboundarycondition_advection_diffusion_2d_gpu - endinterface - - interface - subroutine setgradientboundarycondition_advection_diffusion_2d_gpu(extboundary,boundary,sideinfo,N,nel,nvar) & - bind(c,name="setgradientboundarycondition_advection_diffusion_2d_gpu") - use iso_c_binding - type(c_ptr),value :: extboundary,boundary,sideinfo - integer(c_int),value :: N,nel,nvar - endsubroutine setgradientboundarycondition_advection_diffusion_2d_gpu - endinterface - interface subroutine fluxmethod_advection_diffusion_2d_gpu(solution,solutiongradient,flux,u,v,nu,N,nel,nvar) & bind(c,name="fluxmethod_advection_diffusion_2d_gpu") @@ -82,29 +62,6 @@ subroutine boundaryflux_advection_diffusion_2d_gpu(fb,fextb,dfavg,nhat,nscale,fl contains - subroutine setboundarycondition_advection_diffusion_2d(this) - !! Boundary conditions are set to periodic boundary conditions - implicit none - class(advection_diffusion_2d),intent(inout) :: this - - call setboundarycondition_advection_diffusion_2d_gpu(this%solution%extboundary_gpu, & - this%solution%boundary_gpu,this%mesh%sideInfo_gpu,this%solution%interp%N, & - this%solution%nelem,this%solution%nvar) - - endsubroutine setboundarycondition_advection_diffusion_2d - - subroutine setgradientboundarycondition_advection_diffusion_2d(this) - !! Gradient boundary conditions are set to periodic boundary conditions - implicit none - class(advection_diffusion_2d),intent(inout) :: this - - call setgradientboundarycondition_advection_diffusion_2d_gpu( & - this%solutiongradient%extboundary_gpu, & - this%solutiongradient%boundary_gpu,this%mesh%sideInfo_gpu, & - this%solution%interp%N,this%solution%nelem,this%solution%nvar) - - endsubroutine setgradientboundarycondition_advection_diffusion_2d - subroutine fluxmethod_advection_diffusion_2d(this) implicit none class(advection_diffusion_2d),intent(inout) :: this diff --git a/src/gpu/SELF_advection_diffusion_3d.f90 b/src/gpu/SELF_advection_diffusion_3d.f90 index c470a12cd..51bd89dbf 100644 --- a/src/gpu/SELF_advection_diffusion_3d.f90 +++ b/src/gpu/SELF_advection_diffusion_3d.f90 @@ -32,31 +32,11 @@ module self_advection_diffusion_3d type,extends(advection_diffusion_3d_t) :: advection_diffusion_3d contains - procedure :: setboundarycondition => setboundarycondition_advection_diffusion_3d - procedure :: setgradientboundarycondition => setgradientboundarycondition_advection_diffusion_3d procedure :: boundaryflux => boundaryflux_advection_diffusion_3d procedure :: fluxmethod => fluxmethod_advection_diffusion_3d endtype advection_diffusion_3d - interface - subroutine setboundarycondition_advection_diffusion_3d_gpu(extboundary,boundary,sideinfo,N,nel,nvar) & - bind(c,name="setboundarycondition_advection_diffusion_3d_gpu") - use iso_c_binding - type(c_ptr),value :: extboundary,boundary,sideinfo - integer(c_int),value :: N,nel,nvar - endsubroutine setboundarycondition_advection_diffusion_3d_gpu - endinterface - - interface - subroutine setgradientboundarycondition_advection_diffusion_3d_gpu(extboundary,boundary,sideinfo,N,nel,nvar) & - bind(c,name="setgradientboundarycondition_advection_diffusion_3d_gpu") - use iso_c_binding - type(c_ptr),value :: extboundary,boundary,sideinfo - integer(c_int),value :: N,nel,nvar - endsubroutine setgradientboundarycondition_advection_diffusion_3d_gpu - endinterface - interface subroutine fluxmethod_advection_diffusion_3d_gpu(solution,solutiongradient,flux,u,v,w,nu,N,nel,nvar) & bind(c,name="fluxmethod_advection_diffusion_3d_gpu") @@ -81,29 +61,6 @@ subroutine boundaryflux_advection_diffusion_3d_gpu(fb,fextb,dfavg,nhat,nscale,fl contains - subroutine setboundarycondition_advection_diffusion_3d(this) - !! Boundary conditions are set to periodic boundary conditions - implicit none - class(advection_diffusion_3d),intent(inout) :: this - - call setboundarycondition_advection_diffusion_3d_gpu(this%solution%extboundary_gpu, & - this%solution%boundary_gpu,this%mesh%sideInfo_gpu,this%solution%interp%N, & - this%solution%nelem,this%solution%nvar) - - endsubroutine setboundarycondition_advection_diffusion_3d - - subroutine setgradientboundarycondition_advection_diffusion_3d(this) - !! Gradient boundary conditions are set to periodic boundary conditions - implicit none - class(advection_diffusion_3d),intent(inout) :: this - - call setgradientboundarycondition_advection_diffusion_3d_gpu( & - this%solutiongradient%extboundary_gpu, & - this%solutiongradient%boundary_gpu,this%mesh%sideInfo_gpu, & - this%solution%interp%N,this%solution%nelem,this%solution%nvar) - - endsubroutine setgradientboundarycondition_advection_diffusion_3d - subroutine fluxmethod_advection_diffusion_3d(this) implicit none class(advection_diffusion_3d),intent(inout) :: this diff --git a/test/CMakeLists.txt b/test/CMakeLists.txt index daecac1fe..19a7f72f2 100644 --- a/test/CMakeLists.txt +++ b/test/CMakeLists.txt @@ -128,10 +128,12 @@ add_fortran_tests ( "advection_diffusion_2d_rk3.f90" "advection_diffusion_2d_rk3_pickup.f90" "advection_diffusion_2d_rk4.f90" + "advection_diffusion_2d_nostress.f90" "advection_diffusion_3d_euler.f90" "advection_diffusion_3d_rk2.f90" "advection_diffusion_3d_rk3.f90" "advection_diffusion_3d_rk4.f90" + "advection_diffusion_3d_nostress.f90" "linear_shallow_water_2d_constant.f90" "linear_shallow_water_2d_nonormalflow.f90" "linear_shallow_water_2d_radiation.f90" diff --git a/test/advection_diffusion_2d_nostress.f90 b/test/advection_diffusion_2d_nostress.f90 new file mode 100644 index 000000000..33b58e155 --- /dev/null +++ b/test/advection_diffusion_2d_nostress.f90 @@ -0,0 +1,116 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +program advection_diffusion_2d_nostress +!! Test for the no-normal-flow (mirror) + no-stress (gradient reflection) +!! boundary conditions on the 2D advection-diffusion model. +!! +!! A Gaussian pulse centered in the domain diffuses outward with +!! insulating (zero-flux) walls. The entropy (L2 norm) must decay +!! monotonically, verifying stability of the BR1 no-stress BC. + + use self_data + use self_advection_diffusion_2d + + implicit none + character(SELF_INTEGRATOR_LENGTH),parameter :: integrator = 'rk3' + integer,parameter :: controlDegree = 7 + integer,parameter :: targetDegree = 16 + real(prec),parameter :: nu = 0.005_prec ! diffusivity + real(prec),parameter :: dt = 1.0_prec*10.0_prec**(-4) ! time-step size + real(prec),parameter :: endtime = 0.2_prec + real(prec),parameter :: iointerval = 0.1_prec + real(prec) :: e0,ef ! Initial and final entropy + type(advection_diffusion_2d) :: modelobj + type(Lagrange),target :: interp + type(Mesh2D),target :: mesh + type(SEMQuad),target :: geometry + character(LEN=255) :: WORKSPACE + + ! Create a uniform block mesh + call get_environment_variable("WORKSPACE",WORKSPACE) + call mesh%Read_HOPr(trim(WORKSPACE)//"/share/mesh/Block2D/Block2D_mesh.h5") + + ! Tag all physical boundaries as no-normal-flow walls. + ! This activates both the hyperbolic mirror BC and the + ! parabolic no-stress BC registered by the model. + call mesh%ResetBoundaryConditionType(SELF_BC_NONORMALFLOW) + + ! Create an interpolant + call interp%Init(N=controlDegree, & + controlNodeType=GAUSS, & + M=targetDegree, & + targetNodeType=UNIFORM) + + ! Generate geometry (metric terms) from the mesh elements + call geometry%Init(interp,mesh%nElem) + call geometry%GenerateFromMesh(mesh) + + ! Initialize the model + call modelobj%Init(mesh,geometry) + modelobj%gradient_enabled = .true. + + ! Pure diffusion: zero advection velocity + modelobj%u = 0.0_prec + modelobj%v = 0.0_prec + ! Set the diffusivity + modelobj%nu = nu + + ! Set the initial condition: Gaussian pulse centered in the domain + call modelobj%solution%SetEquation(1,'f = exp( -( (x-0.5)^2 + (y-0.5)^2 )/0.005 )') + call modelobj%solution%SetInteriorFromEquation(geometry,0.0_prec) + + print*,"min, max (interior)", & + minval(modelobj%solution%interior), & + maxval(modelobj%solution%interior) + + call modelobj%CalculateEntropy() + call modelobj%ReportEntropy() + e0 = modelobj%entropy + + ! Set the model's time integration method + call modelobj%SetTimeIntegrator(integrator) + + ! Forward step the model + call modelobj%ForwardStep(endtime,dt,iointerval) + + print*,"min, max (interior)", & + minval(modelobj%solution%interior), & + maxval(modelobj%solution%interior) + ef = modelobj%entropy + + if(ef > e0) then + print*,"Error: Final entropy greater than initial entropy! ",e0,ef + stop 1 + endif + + ! Clean up + call modelobj%free() + call mesh%free() + call geometry%free() + call interp%free() + +endprogram advection_diffusion_2d_nostress diff --git a/test/advection_diffusion_3d_nostress.f90 b/test/advection_diffusion_3d_nostress.f90 new file mode 100644 index 000000000..10c93addd --- /dev/null +++ b/test/advection_diffusion_3d_nostress.f90 @@ -0,0 +1,117 @@ +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! +! +! Maintainers : support@fluidnumerics.com +! Official Repository : https://github.com/FluidNumerics/self/ +! +! Copyright © 2024 Fluid Numerics LLC +! +! Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met: +! +! 1. Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer. +! +! 2. Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in +! the documentation and/or other materials provided with the distribution. +! +! 3. Neither the name of the copyright holder nor the names of its contributors may be used to endorse or promote products derived from +! this software without specific prior written permission. +! +! THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT +! LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT +! HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT +! LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY +! THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF +! THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. +! +! //////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////// ! + +program advection_diffusion_3d_nostress +!! Test for the no-normal-flow (mirror) + no-stress (gradient reflection) +!! boundary conditions on the 3D advection-diffusion model. +!! +!! A Gaussian pulse centered in the domain diffuses outward with +!! insulating (zero-flux) walls. The entropy (L2 norm) must decay +!! monotonically, verifying stability of the BR1 no-stress BC. + + use self_data + use self_advection_diffusion_3d + + implicit none + character(SELF_INTEGRATOR_LENGTH),parameter :: integrator = 'rk3' + integer,parameter :: controlDegree = 7 + integer,parameter :: targetDegree = 16 + real(prec),parameter :: nu = 0.001_prec ! diffusivity + real(prec),parameter :: dt = 1.0_prec*10.0_prec**(-4) ! time-step size + real(prec),parameter :: endtime = 0.01_prec + real(prec),parameter :: iointerval = 0.01_prec + type(advection_diffusion_3d) :: modelobj + type(Lagrange),target :: interp + type(Mesh3D),target :: mesh + type(SEMHex),target :: geometry + character(LEN=255) :: WORKSPACE + real(prec) :: e0,ef + + ! Create a uniform block mesh + call get_environment_variable("WORKSPACE",WORKSPACE) + call mesh%Read_HOPr(trim(WORKSPACE)//"/share/mesh/Block3D/Block3D_mesh.h5") + + ! Tag all physical boundaries as no-normal-flow walls. + ! This activates both the hyperbolic mirror BC and the + ! parabolic no-stress BC registered by the model. + call mesh%ResetBoundaryConditionType(SELF_BC_NONORMALFLOW) + + ! Create an interpolant + call interp%Init(N=controlDegree, & + controlNodeType=GAUSS, & + M=targetDegree, & + targetNodeType=UNIFORM) + + ! Generate geometry (metric terms) from the mesh elements + call geometry%Init(interp,mesh%nElem) + call geometry%GenerateFromMesh(mesh) + + ! Initialize the model + call modelobj%Init(mesh,geometry) + modelobj%gradient_enabled = .true. + + ! Pure diffusion: zero advection velocity + modelobj%u = 0.0_prec + modelobj%v = 0.0_prec + modelobj%w = 0.0_prec + ! Set the diffusivity + modelobj%nu = nu + + ! Set the initial condition: Gaussian pulse centered in the domain + call modelobj%solution%SetEquation(1,'f = exp( -( (x-0.5)^2 + (y-0.5)^2 + (z-0.5)^2 )/0.005 )') + call modelobj%solution%SetInteriorFromEquation(geometry,0.0_prec) + + print*,"min, max (interior)", & + minval(modelobj%solution%interior), & + maxval(modelobj%solution%interior) + + call modelobj%CalculateEntropy() + call modelobj%ReportEntropy() + e0 = modelobj%entropy + + ! Set the model's time integration method + call modelobj%SetTimeIntegrator(integrator) + + ! Forward step the model + call modelobj%ForwardStep(endtime,dt,iointerval) + + print*,"min, max (interior)", & + minval(modelobj%solution%interior), & + maxval(modelobj%solution%interior) + ef = modelobj%entropy + + if(ef > e0) then + print*,"Error: Final entropy greater than initial entropy! ",e0,ef + stop 1 + endif + + ! Clean up + call modelobj%free() + call mesh%free() + call geometry%free() + call interp%free() + +endprogram advection_diffusion_3d_nostress diff --git a/test/ec_advection_2d_entropy_conservation.f90 b/test/ec_advection_2d_entropy_conservation.f90 index 62f06d611..c70e7e553 100644 --- a/test/ec_advection_2d_entropy_conservation.f90 +++ b/test/ec_advection_2d_entropy_conservation.f90 @@ -31,7 +31,7 @@ program ec_advection_2d_entropy_conservation !! Setup: !! - Advection velocity: (u, 0) — purely horizontal !! - Mesh: uniform structured, all boundaries set to SELF_BC_NONORMALFLOW - !! - BC override: hbc2d_NoNormalFlow returns sR = sL (mirror) + !! - BC: ECAdvection2D_t registers a no-normal-flow BC that mirrors sR = sL !! !! With a purely horizontal advection velocity: !! - North/South faces have un = v*ny = 0 (no normal flux at all) @@ -81,7 +81,7 @@ program ec_advection_2d_entropy_conservation targetNodeType=UNIFORM) ! Uniform 5x5 structured mesh on [0,1]x[0,1] with no-normal-flow on all sides. - ! ECAdvection2D_t overrides hbc2d_NoNormalFlow to mirror the interior state, + ! ECAdvection2D_t registers a no-normal-flow BC that mirrors the interior state, ! so sR = sL at every domain face — no upwind dissipation. bcids(1:4) = [SELF_BC_NONORMALFLOW,SELF_BC_NONORMALFLOW, & SELF_BC_NONORMALFLOW,SELF_BC_NONORMALFLOW] diff --git a/test/ec_advection_2d_rk3.f90 b/test/ec_advection_2d_rk3.f90 index 9636d5c02..7efc5639a 100644 --- a/test/ec_advection_2d_rk3.f90 +++ b/test/ec_advection_2d_rk3.f90 @@ -60,7 +60,7 @@ program ec_advection_2d_rk3 targetNodeType=UNIFORM) ! Structured mesh with no-normal-flow BCs on all sides. - ! ECAdvection2D_t overrides hbc2d_NoNormalFlow to mirror (sR=sL), + ! ECAdvection2D_t registers a no-normal-flow BC that mirrors (sR=sL), ! so the upwind Riemann flux has zero dissipation at domain faces. bcids(1:4) = [SELF_BC_NONORMALFLOW,SELF_BC_NONORMALFLOW, & SELF_BC_NONORMALFLOW,SELF_BC_NONORMALFLOW] diff --git a/test/ec_advection_3d_entropy_conservation.f90 b/test/ec_advection_3d_entropy_conservation.f90 index a78eb460d..54af2710a 100644 --- a/test/ec_advection_3d_entropy_conservation.f90 +++ b/test/ec_advection_3d_entropy_conservation.f90 @@ -31,7 +31,7 @@ program ec_advection_3d_entropy_conservation !! Setup: !! - Advection velocity: (u, 0, 0) — purely x-directed !! - Mesh: uniform structured 3x3x3, all boundaries SELF_BC_NONORMALFLOW - !! - BC override: hbc3d_NoNormalFlow returns sR = sL (mirror) + !! - BC: ECAdvection3D_t registers a no-normal-flow BC that mirrors sR = sL !! !! With purely x-directed advection: !! - South/North (y-faces) and Bottom/Top (z-faces) have un = 0