diff --git a/examples/fluids/navierstokes.h b/examples/fluids/navierstokes.h index de240255d9..0c9e575873 100644 --- a/examples/fluids/navierstokes.h +++ b/examples/fluids/navierstokes.h @@ -140,6 +140,15 @@ struct CeedData_private { qf_apply_outflow_jacobian, qf_apply_freestream, qf_apply_freestream_jacobian; }; +// Distance functions +struct Distance_private { + DM dm; + SNES snesDist; + CeedVector q_data; + CeedOperator op_distance_function; + CeedQFunction qf_distance_function; +}; + typedef struct { DM dm; PetscSF sf; // For communicating child data to parents diff --git a/examples/fluids/qfunctions/wall_dist_func.h b/examples/fluids/qfunctions/wall_dist_func.h new file mode 100644 index 0000000000..f3dee8f4dd --- /dev/null +++ b/examples/fluids/qfunctions/wall_dist_func.h @@ -0,0 +1,115 @@ +// Copyright (c) 2017-2023, Lawrence Livermore National Security, LLC and other CEED contributors. +// All Rights Reserved. See the top-level LICENSE and NOTICE files for details. +// +// SPDX-License-Identifier: BSD-2-Clause +// +// This file is part of CEED: http://github.com/ceed + +/// @file + +#ifndef wall_dist_func_h +#define wall_dist_func_h + +#include +#include + +#include "newtonian_state.h" +#include "utils.h" + +CEED_QFUNCTION(SetupDistanceFunctionGeo)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { + // Inputs + const CeedScalar(*J)[3][CEED_Q_VLA] = (const CeedScalar(*)[3][CEED_Q_VLA])in[1]; + const CeedScalar(*w) = in[2]; // Note: *X = in[0] + // Outputs + CeedScalar(*qd) = out[0]; + + const CeedInt dim = 3; + // Quadrature Point Loop + CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { + // Setup + CeedScalar A[3][3]; + for (CeedInt j = 0; j < dim; j++) { + for (CeedInt k = 0; k < dim; k++) { + // Equivalent code with no mod operations: + // A[k][j] = J[k+1][j+1]*J[k+2][j+2] - J[k+1][j+2]*J[k+2][j+1] + A[k][j] = J[(k + 1) % dim][(j + 1) % dim][i] * J[(k + 2) % dim][(j + 2) % dim][i] - + J[(k + 1) % dim][(j + 2) % dim][i] * J[(k + 2) % dim][(j + 1) % dim][i]; + } + } + const CeedScalar detJ = J[0][0][i] * A[0][0] + J[0][1][i] * A[0][1] + J[0][2][i] * A[0][2]; + + const CeedScalar qw = w[i] / detJ; + qd[i + Q * 0] = w[i] * detJ; + qd[i + Q * 1] = qw * (A[0][0] * A[0][0] + A[0][1] * A[0][1] + A[0][2] * A[0][2]); + qd[i + Q * 2] = qw * (A[0][0] * A[1][0] + A[0][1] * A[1][1] + A[0][2] * A[1][2]); + qd[i + Q * 3] = qw * (A[0][0] * A[2][0] + A[0][1] * A[2][1] + A[0][2] * A[2][2]); + qd[i + Q * 4] = qw * (A[1][0] * A[1][0] + A[1][1] * A[1][1] + A[1][2] * A[1][2]); + qd[i + Q * 5] = qw * (A[1][0] * A[2][0] + A[1][1] * A[2][1] + A[1][2] * A[2][2]); + qd[i + Q * 6] = qw * (A[2][0] * A[2][0] + A[2][1] * A[2][1] + A[2][2] * A[2][2]); + } // End of Quadrature Point Loop + return 0; +} + +CEED_QFUNCTION(SetupDistanceFunctionRhs)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { +#ifndef M_PI +#define M_PI 3.14159265358979323846 +#endif + const CeedScalar *x = in[0], *w = in[1]; + CeedScalar *true_soln = out[0], *rhs = out[1]; + + // Quadrature Point Loop + CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { + const CeedScalar c[3] = {0, 1., 2.}; + const CeedScalar k[3] = {1., 2., 3.}; + + true_soln[i] = sin(M_PI * (c[0] + k[0] * x[i + Q * 0])) * sin(M_PI * (c[1] + k[1] * x[i + Q * 1])) * sin(M_PI * (c[2] + k[2] * x[i + Q * 2])); + + rhs[i] = w[i + Q * 0] * M_PI * M_PI * (k[0] * k[0] + k[1] * k[1] + k[2] * k[2]) * true_soln[i]; + } // End of Quadrature Point Loop + + return 0; +} + +// Qfunction for the mass matix +CEED_QFUNCTION(Mass_N)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { + // Inputs + const CeedInt N = 5; + const CeedScalar(*u)[CEED_Q_VLA] = (const CeedScalar(*)[CEED_Q_VLA])in[0]; + const CeedScalar(*q_data) = in[1]; + // const CeedVector ones_vec = 1.0 ; + + // Outputs + CeedScalar(*v)[CEED_Q_VLA] = (CeedScalar(*)[CEED_Q_VLA])out[0]; + + CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { + CeedPragmaSIMD for (CeedInt j = 0; j < N; j++) { v[j][i] = q_data[i] * u[j][i]; } + } + return 0; +} + +CEED_QFUNCTION(DistanceFunction)(void *ctx, CeedInt Q, const CeedScalar *const *in, CeedScalar *const *out) { + // inputs + const CeedScalar *phig = in[0]; + const CeedScalar(*q_data) = in[1]; + // Outputs + CeedScalar *vg = out[0]; + + // Quadrature point loop + CeedPragmaSIMD for (CeedInt i = 0; i < Q; i++) { + // Read spatial derivatives of phi + const CeedScalar dphi[3] = {phig[i + Q * 0], phig[i + Q * 1], phig[i + Q * 2]}; + // Read q_data (dXdxdXdx_T symmetric matrix) + const CeedScalar dXdxdXdx_T[3][3] = { + {q_data[i + 1 * Q], q_data[i + 2 * Q], q_data[i + 3 * Q]}, + {q_data[i + 2 * Q], q_data[i + 4 * Q], q_data[i + 5 * Q]}, + {q_data[i + 3 * Q], q_data[i + 5 * Q], q_data[i + 6 * Q]} + }; + // CeedPragmaSIMD for (CeedInt j = 0; j < 3; j++) { v[j][i] = q_data[i] * u[j][i]; } + for (int j = 0; j < 3; j++) { // j = direction of vg + vg[i + j * Q] = (dphi[0] * dXdxdXdx_T[0][j] + dphi[1] * dXdxdXdx_T[1][j] + dphi[2] * dXdxdXdx_T[2][j]); + } + } + return 0; +} + +#endif // wall_dist_func_h \ No newline at end of file diff --git a/examples/fluids/src/setupdm.c b/examples/fluids/src/setupdm.c index a6668e9170..b610a5102a 100644 --- a/examples/fluids/src/setupdm.c +++ b/examples/fluids/src/setupdm.c @@ -42,12 +42,15 @@ PetscErrorCode SetUpDM(DM dm, ProblemData *problem, PetscInt degree, SimpleBC bc PetscFE fe; PetscInt num_comp_q = 5; DMLabel label; + DM dmDist; PetscCall(PetscFECreateLagrange(PETSC_COMM_SELF, problem->dim, num_comp_q, PETSC_FALSE, degree, PETSC_DECIDE, &fe)); PetscCall(PetscObjectSetName((PetscObject)fe, "Q")); PetscCall(DMAddField(dm, NULL, (PetscObject)fe)); PetscCall(DMCreateDS(dm)); PetscCall(DMGetLabel(dm, "Face Sets", &label)); PetscCall(DMPlexLabelComplete(dm, label)); + PetscCall(DMClone(dm, &dmDist)); + PetscCall(DMAddField(dmDist, NULL, (PetscObject)fe)); // Set wall BCs if (bc->num_wall > 0) { PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "wall", label, bc->num_wall, bc->walls, 0, bc->num_comps, bc->wall_comps, @@ -71,6 +74,11 @@ PetscErrorCode SetUpDM(DM dm, ProblemData *problem, PetscInt degree, SimpleBC bc PetscCall(DMAddBoundary(dm, DM_BC_ESSENTIAL, "slipz", label, bc->num_slip[2], bc->slips[2], 0, 1, comps, (void (*)(void))NULL, NULL, problem->bc_ctx, NULL)); } + // Set dirichlet boundary condition on cylinder walls + if (bc->num_wall > 0) { + PetscCall(DMAddBoundary(dmDist, DM_BC_ESSENTIAL, "cylinder_wall", label, bc->num_wall, bc->walls, 0, 1, bc->wall_comps, + (void (*)(void))problem->bc, NULL, problem->bc_ctx, NULL)); + } { PetscBool use_strongstg = PETSC_FALSE; PetscCall(PetscOptionsGetBool(NULL, NULL, "-stg_strong", &use_strongstg, NULL)); diff --git a/examples/fluids/src/wall_dist_func.c b/examples/fluids/src/wall_dist_func.c new file mode 100644 index 0000000000..e5382200ca --- /dev/null +++ b/examples/fluids/src/wall_dist_func.c @@ -0,0 +1,161 @@ +// Copyright (c) 2017-2023, Lawrence Livermore National Security, LLC and other CEED contributors. +// All Rights Reserved. See the top-level LICENSE and NOTICE files for details. +// +// SPDX-License-Identifier: BSD-2-Clause +// +// This file is part of CEED: http://github.com/ceed + +/// @file +/// General wall distance functions for Navier-Stokes example using PETSc +/// We do this by solving the Poisson equation ∇^{2} φ = -1 with weak form ∫ ∇v ⋅ ∇ φ - v 1 + +#include "../qfunctions/wall_dist_func.h" + +#include "../navierstokes.h" +#include "../qfunctions/newtonian_state.h" + +// General distance functions +static PetscErrorCode Distance_Function_NS(DM dm, User user) { + MPI_Comm comm; + PetscScalar *r; + DM dmDist; + PetscFE fe; + PetscInt xl_size, l_size, g_size, dim = 3; + SNES snesDist; + Mat A; + Vec X, Y, q_data, X_loc, rhs, rhs_loc; + VecType vec_type; + PetscMemType mem_type; + Ceed ceed; + CeedInt num_elem = 10, P = 3, Q = P + 1; + CeedInt num_nodes_x = num_elem + 1, num_nodes_phi = num_elem * (P - 1) + 1; + CeedInt ind_x[num_elem * 2], ind_phi[num_elem * P]; + CeedInt strides_qd[3] = {1, Q, Q}; + CeedVector rhs_ceed, X_ceed, Y_ceed, q_data_ceed, ones_vec; + CeedOperator op_distance_function, op_mass; + CeedQFunction qf_distance_function, qf_mass; + CeedBasis basis_x, basis_phi; + CeedElemRestriction elem_restr_x, elem_restr_phi, elem_restr_qd; + + PetscFunctionBeginUser; + comm = PETSC_COMM_WORLD; + PetscCall(PetscFECreateDefault(PETSC_COMM_SELF, dim, 1, PETSC_FALSE, NULL, PETSC_DETERMINE, &fe)); + PetscBool distance_snes_monitor = PETSC_FALSE; + PetscCall(PetscOptionsHasName(NULL, NULL, "-distance_snes_monitor", &distance_snes_monitor)); + PetscCall(SNESCreate(PETSC_COMM_WORLD, &snesDist)); + PetscObjectSetOptionsPrefix((PetscObject)snesDist, "distance_"); + PetscCall(DMClone(dm, &dmDist)); + PetscCall(DMGetVecType(dm, &vec_type)); + PetscCall(DMSetVecType(dmDist, vec_type)); + PetscCall(DMAddField(dmDist, NULL, (PetscObject)fe)); + + // Create Vectors + PetscCall(DMCreateGlobalVector(dmDist, &X)); + PetscCall(DMCreateGlobalVector(dmDist, &Y)); + PetscCall(VecGetLocalSize(X, &l_size)); + PetscCall(VecGetSize(X, &g_size)); + PetscCall(VecGetSize(q_data, &g_size)); + PetscCall(DMCreateLocalVector(dmDist, &X_loc)); + PetscCall(VecGetSize(X_loc, &xl_size)); + PetscCall(VecDuplicate(X, &rhs)); + + // Create Matshell Operator + PetscCall(MatCreateShell(comm, l_size, l_size, g_size, g_size, op_distance_function, &A)); + PetscCall(MatShellSetOperation(A, MATOP_MULT, (void (*)(void))MatMult_Ceed)); + PetscCall(MatShellSetVecType(A, vec_type)); + + // Setup libCEED vector + PetscCall(VecP2C(X, &mem_type, X_ceed)); + PetscCall(VecP2C(Y, &mem_type, Y_ceed)); + PetscCall(VecP2C(q_data, &mem_type, q_data_ceed)); + + // Create RHS vector + PetscCall(VecDuplicate(X_loc, &rhs_loc)); + PetscCall(VecZeroEntries(rhs_loc)); + PetscCall(VecGetArrayAndMemType(rhs_loc, &r, &mem_type)); + CeedVectorCreate(ceed, xl_size, &rhs_ceed); + CeedVectorSetArray(rhs_ceed, MemTypeP2C(mem_type), CEED_USE_POINTER, r); + + // Gather RHS + CeedVectorTakeArray(rhs_ceed, MemTypeP2C(mem_type), NULL); + PetscCall(VecRestoreArrayAndMemType(rhs_loc, &r)); + PetscCall(VecZeroEntries(rhs)); + PetscCall(DMLocalToGlobal(dmDist, rhs_loc, ADD_VALUES, rhs)); + CeedVectorDestroy(&rhs_ceed); + + // Create Element Restriction + CeedElemRestrictionCreate(ceed, num_elem, 2, 1, 1, num_nodes_x, CEED_MEM_HOST, CEED_USE_POINTER, ind_x, &elem_restr_x); + CeedElemRestrictionCreate(ceed, num_elem, 2, 1, 1, num_nodes_phi, CEED_MEM_HOST, CEED_USE_POINTER, ind_phi, &elem_restr_phi); + CeedElemRestrictionCreateStrided(ceed, num_elem, Q, 1, Q * num_elem, strides_qd, &elem_restr_qd); + + // Create Basis + CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, P, Q, CEED_GAUSS, &basis_x); + CeedBasisCreateTensorH1Lagrange(ceed, dim, 1, P, Q, CEED_GAUSS, &basis_phi); + + // Create and Add QFunction fields + CeedQFunctionCreateInterior(ceed, dim, DistanceFunction, DistanceFunction_loc, &qf_distance_function); + // CeedQFunctionCreateInterior(ceed, dim, DistanceFunction - Mass_N, DistanceFunction_loc, &qf_distance_function); + CeedQFunctionCreateInterior(ceed, dim, Mass_N, Mass_N_loc, &qf_mass); + + CeedQFunctionAddInput(qf_distance_function, "v", dim, CEED_EVAL_INTERP); + CeedQFunctionAddInput(qf_distance_function, "dphi", dim, CEED_EVAL_GRAD); + CeedQFunctionAddInput(qf_distance_function, "q_data", dim * (dim + 1) / 2, CEED_EVAL_NONE); + CeedQFunctionAddOutput(qf_distance_function, "dv", dim, CEED_EVAL_GRAD); + + CeedQFunctionAddInput(qf_mass, "rho", 1, CEED_EVAL_NONE); + CeedQFunctionAddInput(qf_mass, "phi", 1, CEED_EVAL_INTERP); + CeedQFunctionAddOutput(qf_mass, "v", 1, CEED_EVAL_INTERP); + + // Create Operator + CeedOperatorCreate(ceed, qf_distance_function, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_distance_function); + CeedOperatorCreate(ceed, qf_mass, CEED_QFUNCTION_NONE, CEED_QFUNCTION_NONE, &op_mass); + + // Operator set + CeedOperatorSetField(op_distance_function, "phi", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_distance_function, "v", elem_restr_x, basis_x, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_distance_function, "output", elem_restr_phi, CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE); + + CeedOperatorSetField(op_mass, "rho", elem_restr_qd, CEED_BASIS_COLLOCATED, q_data_ceed); + CeedOperatorSetField(op_mass, "phi", elem_restr_phi, basis_phi, CEED_VECTOR_ACTIVE); + CeedOperatorSetField(op_mass, "v", elem_restr_phi, basis_phi, CEED_VECTOR_ACTIVE); + + CeedVectorCreate(ceed, num_nodes_phi, &X_ceed); + CeedVectorSetValue(X_ceed, 0.0); + CeedVectorCreate(ceed, num_nodes_phi, &Y_ceed); + CeedVectorSetValue(Y_ceed, 0.0); + CeedVectorCreate(ceed, num_nodes_phi, &ones_vec); + CeedVectorSetValue(ones_vec, 1.0); // vector of ones to multiply mass matrix with + + // Apply Operator + CeedOperatorApply(op_distance_function, X_ceed, Y_ceed, CEED_REQUEST_IMMEDIATE); + CeedOperatorApply(op_mass, X_ceed, ones_vec, CEED_REQUEST_IMMEDIATE); + + // Restore PETSc vectors + PetscCall(VecC2P(X_ceed, mem_type, X)); + PetscCall(VecC2P(Y_ceed, mem_type, Y)); + + // Set up SNES + PetscCall(SNESCreate(PETSC_COMM_WORLD, &snesDist)); + + // Solve + PetscCall(SNESSolve(snesDist, rhs, X)); + PetscCall(SNESGetSolution(snesDist, &X)); + + // Clean up + PetscCall(VecDestroy(&rhs)); + PetscCall(VecDestroy(&rhs_loc)); + PetscCall(SNESDestroy(&snesDist)); + PetscCall(VecDestroy(&X)); + PetscCall(VecDestroy(&X_loc)); + PetscCall(VecDestroy(&q_data)); + PetscCall(MatDestroy(&A)); + CeedVectorDestroy(&ones_vec); + CeedElemRestrictionDestroy(&elem_restr_x); + CeedElemRestrictionDestroy(&elem_restr_phi); + CeedElemRestrictionDestroy(&elem_restr_qd); + CeedQFunctionDestroy(&qf_distance_function); + CeedOperatorDestroy(&op_distance_function); + CeedOperatorDestroy(&op_mass); + CeedDestroy(&ceed); + PetscFunctionReturn(0); +}