static char help[] = "Discontinuous Galerkin advection example.\n\n\
This example demonstrates DG method implementation in PETSc:\n\
  - Creating discontinuous finite elements with PetscFECreateBrokenElement()\n\
  - Setting up Riemann solvers for interface fluxes\n\
  - Implementing DG weak form for scalar advection\n\
  - Time integration with SSPRK methods\n\n\
Run with: ./ex3 -dm_plex_dim 2 -dm_plex_box_faces 10,10 -ts_max_steps 100\n\n";

/*
  Solves the 2D advection equation using Discontinuous Galerkin:
    u_t + div(v*u) = 0
  where v = (v_x, v_y) is a constant velocity field.
  
  Initial condition: Gaussian bump
  Boundary conditions: Periodic or inflow/outflow
*/

#include <petscdmplex.h>
#include <petscds.h>
#include <petscts.h>
#include <petscfe.h>

typedef struct {
  PetscReal velocity[2]; /* Advection velocity */
  PetscReal center[2];   /* Initial bump center */
  PetscReal radius;      /* Initial bump radius */
} AppCtx;

/* 
  Riemann solver for scalar advection using upwind flux
  flux = v*u where u is chosen based on upwind direction
*/
static void RiemannSolver_Advection(PetscInt dim, PetscInt Nf, const PetscReal *qp, const PetscReal *n, const PetscScalar *uL, const PetscScalar *uR, PetscInt numConstants, const PetscScalar constants[], PetscScalar *flux, void *ctx)
{
  AppCtx   *user = (AppCtx *)ctx;
  PetscReal vn   = 0.0;
  PetscInt  d;

  /* Compute velocity dot normal: v·n */
  for (d = 0; d < dim; ++d) vn += user->velocity[d] * n[d];

  /* Upwind flux: choose u based on sign of v·n */
  if (vn > 0.0) flux[0] = vn * uL[0]; /* Flow from left */
  else flux[0] = vn * uR[0];          /* Flow from right */
}

/* 
  Volume integral: time derivative term
  f0 = u_t (coefficient of test function v)
*/
static void f0_advection(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f0[])
{
  f0[0] = u_t[0];
}

/* 
  Volume integral: flux divergence term
  f1 = -v*u (coefficient of grad(v), the test function gradient)
  Weak form: -∫ grad(v)·(velocity*u) dx
*/
static void f1_advection(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt uOff[], const PetscInt uOff_x[], const PetscScalar u[], const PetscScalar u_t[], const PetscScalar u_x[], const PetscInt aOff[], const PetscInt aOff_x[], const PetscScalar a[], const PetscScalar a_t[], const PetscScalar a_x[], PetscReal t, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar f1[])
{
  AppCtx  *user = (AppCtx *)constants;
  PetscInt d;

  for (d = 0; d < dim; ++d) f1[d] = -user->velocity[d] * u[0];
}

/* 
  Inflow boundary condition: set ghost state to zero
*/
static PetscErrorCode BoundaryInflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *uI, PetscScalar *uG, void *ctx)
{
  PetscFunctionBeginUser;
  uG[0] = 0.0; /* Inflow state */
  PetscFunctionReturn(PETSC_SUCCESS);
}

/* 
  Outflow boundary condition: set ghost state equal to interior
*/
static PetscErrorCode BoundaryOutflow(PetscReal time, const PetscReal *c, const PetscReal *n, const PetscScalar *uI, PetscScalar *uG, void *ctx)
{
  PetscFunctionBeginUser;
  uG[0] = uI[0]; /* Outflow: copy interior state */
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*
  Initial condition function: Gaussian bump
  u(x,y,0) = exp(-r^2/radius^2) where r = |x - center|
*/
static PetscErrorCode InitialCondition(PetscInt dim, PetscReal t, const PetscReal x[], PetscInt Nf, PetscScalar u[], void *ctx)
{
  AppCtx   *user = (AppCtx *)ctx;
  PetscReal dx, dy, r2;

  PetscFunctionBeginUser;
  dx   = x[0] - user->center[0];
  dy   = x[1] - user->center[1];
  r2   = dx * dx + dy * dy;
  u[0] = PetscExpReal(-r2 / (user->radius * user->radius));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/*
  Set initial condition on the solution vector
*/
static PetscErrorCode SetInitialCondition(DM dm, Vec U, AppCtx *user)
{
  PetscErrorCode (*func[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
  void *ctx[1];

  PetscFunctionBeginUser;
  func[0] = InitialCondition;
  ctx[0]  = user;

  PetscCall(DMProjectFunction(dm, 0.0, func, ctx, INSERT_ALL_VALUES, U));
  PetscFunctionReturn(PETSC_SUCCESS);
}

/* Monitor function to print solution statistics */
static PetscErrorCode Monitor(TS ts, PetscInt step, PetscReal time, Vec U, void *ctx)
{
  DM        dm;
  PetscReal min, max, norm;

  PetscFunctionBeginUser;
  if (step < 0) PetscFunctionReturn(PETSC_SUCCESS); /* Step -1 is before solving */

  PetscCall(TSGetDM(ts, &dm));
  PetscCall(VecMin(U, NULL, &min));
  PetscCall(VecMax(U, NULL, &max));
  PetscCall(VecNorm(U, NORM_2, &norm));

  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Step %4" PetscInt_FMT " time %8.4f  min %8.4e  max %8.4e  L2norm %8.4e\n", step, (double)time, (double)min, (double)max, (double)norm));
  PetscFunctionReturn(PETSC_SUCCESS);
}

int main(int argc, char **argv)
{
  DM             dm;
  TS             ts;
  Vec            u;
  PetscFE        fe, dgfe;
  PetscDS        ds;
  AppCtx         user;
  PetscInt       dim = 2, order = 1;
  PetscBool      simplex = PETSC_FALSE;
  PetscReal      ftime, dt;
  DMLabel        label;
  const PetscInt inflowids[]  = {1, 2}; /* Left and bottom boundaries */
  const PetscInt outflowids[] = {3, 4}; /* Right and top boundaries */

  PetscFunctionBeginUser;
  PetscCall(PetscInitialize(&argc, &argv, NULL, help));

  /* Set default parameters */
  user.velocity[0] = 1.0; /* Advection velocity in x */
  user.velocity[1] = 0.5; /* Advection velocity in y */
  user.center[0]   = 0.25;
  user.center[1]   = 0.25;
  user.radius      = 0.15;

  PetscOptionsBegin(PETSC_COMM_WORLD, "", "DG Advection Options", "");
  PetscCall(PetscOptionsInt("-order", "Polynomial order of DG elements", "", order, &order, NULL));
  PetscCall(PetscOptionsRealArray("-velocity", "Advection velocity (vx,vy)", "", user.velocity, &dim, NULL));
  PetscCall(PetscOptionsReal("-radius", "Radius of initial Gaussian bump", "", user.radius, &user.radius, NULL));
  PetscOptionsEnd();

  /* Create mesh */
  PetscCall(DMCreate(PETSC_COMM_WORLD, &dm));
  PetscCall(DMSetType(dm, DMPLEX));
  PetscCall(DMSetFromOptions(dm));
  PetscCall(DMGetDimension(dm, &dim));
  PetscCall(DMPlexIsSimplex(dm, &simplex));
  PetscCall(DMSetUp(dm));
  PetscCall(DMViewFromOptions(dm, NULL, "-dm_view"));

  /* Create continuous finite element */
  PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, dim, 1, simplex, order, PETSC_DETERMINE, &fe));
  PetscCall(PetscObjectSetName((PetscObject)fe, "potential"));

  /* Convert to DISCONTINUOUS element - this is the key for DG! */
  PetscCall(PetscFECreateBrokenElement(fe, &dgfe));
  PetscCall(PetscFEDestroy(&fe));
  PetscCall(PetscObjectSetName((PetscObject)dgfe, "DG_potential"));

  /* Add discontinuous field to DM */
  PetscCall(DMAddField(dm, NULL, (PetscObject)dgfe));
  PetscCall(DMCreateDS(dm));
  PetscCall(DMGetDS(dm, &ds));

  /* Set up DG weak form */
  /* Volume integrals */
  PetscCall(PetscDSSetResidual(ds, 0, f0_advection, f1_advection));

  /* Interface fluxes via Riemann solver */
  PetscCall(PetscDSSetRiemannSolver(ds, 0, RiemannSolver_Advection));
  PetscCall(PetscDSSetContext(ds, 0, &user));

  /* Boundary conditions - Riemann BC (natural/weak) */
  PetscCall(DMGetLabel(dm, "marker", &label));
  if (!label) PetscCall(DMCreateLabel(dm, "marker"));
  PetscCall(DMGetLabel(dm, "marker", &label));

  /* Set inflow and outflow boundaries */
  PetscCall(PetscDSAddBoundary(ds, DM_BC_NATURAL_RIEMANN, "inflow", label, PETSC_STATIC_ARRAY_LENGTH(inflowids), inflowids, 0, 0, NULL, (void (*)(void))BoundaryInflow, NULL, &user, NULL));
  PetscCall(PetscDSAddBoundary(ds, DM_BC_NATURAL_RIEMANN, "outflow", label, PETSC_STATIC_ARRAY_LENGTH(outflowids), outflowids, 0, 0, NULL, (void (*)(void))BoundaryOutflow, NULL, &user, NULL));

  /* Create time stepper */
  PetscCall(TSCreate(PETSC_COMM_WORLD, &ts));
  PetscCall(TSSetDM(ts, dm));
  PetscCall(TSSetProblemType(ts, TS_NONLINEAR));
  PetscCall(TSSetType(ts, TSSSP));                  /* Strong Stability Preserving Runge-Kutta */
  PetscCall(TSSetMaxTime(ts, 1.0));                 /* Final time */
  PetscCall(TSSetExactFinalTime(ts, TS_EXACTFINALTIME_MATCHSTEP));

  /* Set up implicit RHS function for DG */
  PetscCall(DMTSSetRHSFunctionLocal(dm, DMPlexTSComputeRHSFunctionFVM, &user));

  /* CFL-based time step */
  PetscReal hmin, maxvel;
  PetscCall(DMPlexGetGeometryFVM(dm, NULL, NULL, &hmin));
  maxvel = PetscSqrtReal(user.velocity[0] * user.velocity[0] + user.velocity[1] * user.velocity[1]);
  dt     = 0.5 * hmin / maxvel; /* CFL condition */
  PetscCall(TSSetTimeStep(ts, dt));

  PetscCall(TSMonitorSet(ts, Monitor, &user, NULL));
  PetscCall(TSSetFromOptions(ts));

  /* Create solution vector and set initial condition */
  PetscCall(DMCreateGlobalVector(dm, &u));
  PetscCall(PetscObjectSetName((PetscObject)u, "solution"));
  PetscCall(SetInitialCondition(dm, u, &user));

  /* Solve */
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nSolving 2D DG advection with velocity (%g, %g)\n", (double)user.velocity[0], (double)user.velocity[1]));
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "DG polynomial order: %" PetscInt_FMT "\n", order));
  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Time step (CFL): %g\n\n", (double)dt));

  PetscCall(TSSolve(ts, u));
  PetscCall(TSGetSolveTime(ts, &ftime));

  PetscCall(PetscPrintf(PETSC_COMM_WORLD, "\nFinal time: %g\n", (double)ftime));

  /* Visualize solution */
  PetscCall(VecViewFromOptions(u, NULL, "-sol_view"));

  /* Clean up */
  PetscCall(VecDestroy(&u));
  PetscCall(TSDestroy(&ts));
  PetscCall(PetscFEDestroy(&dgfe));
  PetscCall(DMDestroy(&dm));
  PetscCall(PetscFinalize());
  return 0;
}

/*TEST

  testset:
    args: -dm_plex_box_faces 8,8 -dm_plex_dim 2 -dm_plex_simplex 0 \
          -order 1 -ts_max_steps 10 -ts_monitor

    test:
      suffix: 2d_quad_p1
      args: -velocity 1.0,0.5

    test:
      suffix: 2d_quad_p2
      args: -order 2 -velocity 1.0,0.0

    test:
      suffix: 2d_tri_p1
      requires: triangle
      args: -dm_plex_simplex 1 -velocity 1.0,1.0

TEST*/