bps.c (288c044332e33f37503f09b6484fec9d0a55fba1) bps.c (a2fa791084b022ffd3d077ae9be76b4b170bb84e)
1// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
2// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
3// reserved. See files LICENSE and NOTICE for details.
4//
5// This file is part of CEED, a collection of benchmarks, miniapps, software
6// libraries and APIs for efficient high-order finite element and spectral
7// element discretizations for exascale applications. For more information and
8// source code availability see http://github.com/ceed.

--- 124 unchanged lines hidden (view full) ---

133struct User_ {
134 MPI_Comm comm;
135 VecScatter ltog; // Scatter for all entries
136 VecScatter ltog0; // Skip Dirichlet values
137 VecScatter gtogD; // global-to-global; only Dirichlet values
138 Vec Xloc, Yloc;
139 CeedVector xceed, yceed;
140 CeedOperator op;
1// Copyright (c) 2017, Lawrence Livermore National Security, LLC. Produced at
2// the Lawrence Livermore National Laboratory. LLNL-CODE-734707. All Rights
3// reserved. See files LICENSE and NOTICE for details.
4//
5// This file is part of CEED, a collection of benchmarks, miniapps, software
6// libraries and APIs for efficient high-order finite element and spectral
7// element discretizations for exascale applications. For more information and
8// source code availability see http://github.com/ceed.

--- 124 unchanged lines hidden (view full) ---

133struct User_ {
134 MPI_Comm comm;
135 VecScatter ltog; // Scatter for all entries
136 VecScatter ltog0; // Skip Dirichlet values
137 VecScatter gtogD; // global-to-global; only Dirichlet values
138 Vec Xloc, Yloc;
139 CeedVector xceed, yceed;
140 CeedOperator op;
141 CeedVector rho;
141 CeedVector qdata;
142 Ceed ceed;
143};
144
145// BP Options
146typedef enum {
147 CEED_BP1 = 0, CEED_BP2 = 1, CEED_BP3 = 2,
148 CEED_BP4 = 3, CEED_BP5 = 4, CEED_BP6 = 5
149} bpType;

--- 242 unchanged lines hidden (view full) ---

392 VecScatter ltog, ltog0, gtogD;
393 User user;
394 Ceed ceed;
395 CeedBasis basisx, basisu;
396 CeedElemRestriction Erestrictx, Erestrictu, Erestrictxi, Erestrictui,
397 Erestrictqdi;
398 CeedQFunction qf_setup, qf_apply, qf_error;
399 CeedOperator op_setup, op_apply, op_error;
142 Ceed ceed;
143};
144
145// BP Options
146typedef enum {
147 CEED_BP1 = 0, CEED_BP2 = 1, CEED_BP3 = 2,
148 CEED_BP4 = 3, CEED_BP5 = 4, CEED_BP6 = 5
149} bpType;

--- 242 unchanged lines hidden (view full) ---

392 VecScatter ltog, ltog0, gtogD;
393 User user;
394 Ceed ceed;
395 CeedBasis basisx, basisu;
396 CeedElemRestriction Erestrictx, Erestrictu, Erestrictxi, Erestrictui,
397 Erestrictqdi;
398 CeedQFunction qf_setup, qf_apply, qf_error;
399 CeedOperator op_setup, op_apply, op_error;
400 CeedVector xcoord, rho, rhsceed, target;
400 CeedVector xcoord, qdata, rhsceed, target;
401 CeedInt P, Q;
402 const CeedInt dim = 3, ncompx = 3;
403 bpType bpChoice;
404
405 ierr = PetscInitialize(&argc, &argv, NULL, help);
406 if (ierr) return ierr;
407 comm = PETSC_COMM_WORLD;
408 ierr = PetscOptionsBegin(comm, NULL, "CEED BPs in PETSc", NULL); CHKERRQ(ierr);

--- 210 unchanged lines hidden (view full) ---

619
620 // Create the Q-function that builds the operator (i.e. computes its
621 // quadrature data) and set its context data
622 CeedQFunctionCreateInterior(ceed, 1, bpOptions[bpChoice].setup,
623 bpOptions[bpChoice].setupfname, &qf_setup);
624 CeedQFunctionAddInput(qf_setup, "x", ncompx, CEED_EVAL_INTERP);
625 CeedQFunctionAddInput(qf_setup, "dx", ncompx*dim, CEED_EVAL_GRAD);
626 CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT);
401 CeedInt P, Q;
402 const CeedInt dim = 3, ncompx = 3;
403 bpType bpChoice;
404
405 ierr = PetscInitialize(&argc, &argv, NULL, help);
406 if (ierr) return ierr;
407 comm = PETSC_COMM_WORLD;
408 ierr = PetscOptionsBegin(comm, NULL, "CEED BPs in PETSc", NULL); CHKERRQ(ierr);

--- 210 unchanged lines hidden (view full) ---

619
620 // Create the Q-function that builds the operator (i.e. computes its
621 // quadrature data) and set its context data
622 CeedQFunctionCreateInterior(ceed, 1, bpOptions[bpChoice].setup,
623 bpOptions[bpChoice].setupfname, &qf_setup);
624 CeedQFunctionAddInput(qf_setup, "x", ncompx, CEED_EVAL_INTERP);
625 CeedQFunctionAddInput(qf_setup, "dx", ncompx*dim, CEED_EVAL_GRAD);
626 CeedQFunctionAddInput(qf_setup, "weight", 1, CEED_EVAL_WEIGHT);
627 CeedQFunctionAddOutput(qf_setup, "rho", bpOptions[bpChoice].qdatasize,
627 CeedQFunctionAddOutput(qf_setup, "qdata", bpOptions[bpChoice].qdatasize,
628 CEED_EVAL_NONE);
629 CeedQFunctionAddOutput(qf_setup, "true_soln", ncompu, CEED_EVAL_NONE);
630 CeedQFunctionAddOutput(qf_setup, "rhs", ncompu, CEED_EVAL_INTERP);
631
632 // Set up PDE operator
633 CeedQFunctionCreateInterior(ceed, 1, bpOptions[bpChoice].apply,
634 bpOptions[bpChoice].applyfname, &qf_apply);
635 // Add inputs and outputs
636 CeedInt gradInScale = bpOptions[bpChoice].inmode==CEED_EVAL_GRAD ? 3 : 1;
637 CeedInt gradOutScale = bpOptions[bpChoice].outmode==CEED_EVAL_GRAD ? 3 : 1;
638 CeedQFunctionAddInput(qf_apply, "u", ncompu*gradInScale,
639 bpOptions[bpChoice].inmode);
628 CEED_EVAL_NONE);
629 CeedQFunctionAddOutput(qf_setup, "true_soln", ncompu, CEED_EVAL_NONE);
630 CeedQFunctionAddOutput(qf_setup, "rhs", ncompu, CEED_EVAL_INTERP);
631
632 // Set up PDE operator
633 CeedQFunctionCreateInterior(ceed, 1, bpOptions[bpChoice].apply,
634 bpOptions[bpChoice].applyfname, &qf_apply);
635 // Add inputs and outputs
636 CeedInt gradInScale = bpOptions[bpChoice].inmode==CEED_EVAL_GRAD ? 3 : 1;
637 CeedInt gradOutScale = bpOptions[bpChoice].outmode==CEED_EVAL_GRAD ? 3 : 1;
638 CeedQFunctionAddInput(qf_apply, "u", ncompu*gradInScale,
639 bpOptions[bpChoice].inmode);
640 CeedQFunctionAddInput(qf_apply, "rho", bpOptions[bpChoice].qdatasize,
640 CeedQFunctionAddInput(qf_apply, "qdata", bpOptions[bpChoice].qdatasize,
641 CEED_EVAL_NONE);
642 CeedQFunctionAddOutput(qf_apply, "v", ncompu*gradOutScale,
643 bpOptions[bpChoice].outmode);
644
645 // Create the error qfunction
646 CeedQFunctionCreateInterior(ceed, 1, bpOptions[bpChoice].error,
647 bpOptions[bpChoice].errorfname, &qf_error);
648 CeedQFunctionAddInput(qf_error, "u", ncompu, CEED_EVAL_INTERP);
649 CeedQFunctionAddInput(qf_error, "true_soln", ncompu, CEED_EVAL_NONE);
650 CeedQFunctionAddOutput(qf_error, "error", ncompu, CEED_EVAL_NONE);
651
652 // Create the persistent vectors that will be needed in setup
653 CeedInt nqpts;
654 CeedBasisGetNumQuadraturePoints(basisu, &nqpts);
641 CEED_EVAL_NONE);
642 CeedQFunctionAddOutput(qf_apply, "v", ncompu*gradOutScale,
643 bpOptions[bpChoice].outmode);
644
645 // Create the error qfunction
646 CeedQFunctionCreateInterior(ceed, 1, bpOptions[bpChoice].error,
647 bpOptions[bpChoice].errorfname, &qf_error);
648 CeedQFunctionAddInput(qf_error, "u", ncompu, CEED_EVAL_INTERP);
649 CeedQFunctionAddInput(qf_error, "true_soln", ncompu, CEED_EVAL_NONE);
650 CeedQFunctionAddOutput(qf_error, "error", ncompu, CEED_EVAL_NONE);
651
652 // Create the persistent vectors that will be needed in setup
653 CeedInt nqpts;
654 CeedBasisGetNumQuadraturePoints(basisu, &nqpts);
655 CeedVectorCreate(ceed, bpOptions[bpChoice].qdatasize*nelem*nqpts, &rho);
655 CeedVectorCreate(ceed, bpOptions[bpChoice].qdatasize*nelem*nqpts, &qdata);
656 CeedVectorCreate(ceed, nelem*nqpts*ncompu, &target);
657 CeedVectorCreate(ceed, lsize*ncompu, &rhsceed);
658
659 // Create the operator that builds the quadrature data for the ceed operator
660 CeedOperatorCreate(ceed, qf_setup, NULL, NULL, &op_setup);
661 CeedOperatorSetField(op_setup, "x", Erestrictx, CEED_NOTRANSPOSE,
662 basisx, CEED_VECTOR_ACTIVE);
663 CeedOperatorSetField(op_setup, "dx", Erestrictx, CEED_NOTRANSPOSE,
664 basisx, CEED_VECTOR_ACTIVE);
665 CeedOperatorSetField(op_setup, "weight", Erestrictxi, CEED_NOTRANSPOSE,
666 basisx, CEED_VECTOR_NONE);
656 CeedVectorCreate(ceed, nelem*nqpts*ncompu, &target);
657 CeedVectorCreate(ceed, lsize*ncompu, &rhsceed);
658
659 // Create the operator that builds the quadrature data for the ceed operator
660 CeedOperatorCreate(ceed, qf_setup, NULL, NULL, &op_setup);
661 CeedOperatorSetField(op_setup, "x", Erestrictx, CEED_NOTRANSPOSE,
662 basisx, CEED_VECTOR_ACTIVE);
663 CeedOperatorSetField(op_setup, "dx", Erestrictx, CEED_NOTRANSPOSE,
664 basisx, CEED_VECTOR_ACTIVE);
665 CeedOperatorSetField(op_setup, "weight", Erestrictxi, CEED_NOTRANSPOSE,
666 basisx, CEED_VECTOR_NONE);
667 CeedOperatorSetField(op_setup, "rho", Erestrictqdi, CEED_NOTRANSPOSE,
667 CeedOperatorSetField(op_setup, "qdata", Erestrictqdi, CEED_NOTRANSPOSE,
668 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
669 CeedOperatorSetField(op_setup, "true_soln", Erestrictui, CEED_NOTRANSPOSE,
670 CEED_BASIS_COLLOCATED, target);
671 CeedOperatorSetField(op_setup, "rhs", Erestrictu, CEED_TRANSPOSE,
672 basisu, rhsceed);
673
674 // Create the mass or diff operator
675 CeedOperatorCreate(ceed, qf_apply, NULL, NULL, &op_apply);
676 CeedOperatorSetField(op_apply, "u", Erestrictu, CEED_TRANSPOSE,
677 basisu, CEED_VECTOR_ACTIVE);
668 CEED_BASIS_COLLOCATED, CEED_VECTOR_ACTIVE);
669 CeedOperatorSetField(op_setup, "true_soln", Erestrictui, CEED_NOTRANSPOSE,
670 CEED_BASIS_COLLOCATED, target);
671 CeedOperatorSetField(op_setup, "rhs", Erestrictu, CEED_TRANSPOSE,
672 basisu, rhsceed);
673
674 // Create the mass or diff operator
675 CeedOperatorCreate(ceed, qf_apply, NULL, NULL, &op_apply);
676 CeedOperatorSetField(op_apply, "u", Erestrictu, CEED_TRANSPOSE,
677 basisu, CEED_VECTOR_ACTIVE);
678 CeedOperatorSetField(op_apply, "rho", Erestrictqdi, CEED_NOTRANSPOSE,
679 CEED_BASIS_COLLOCATED, rho);
678 CeedOperatorSetField(op_apply, "qdata", Erestrictqdi, CEED_NOTRANSPOSE,
679 CEED_BASIS_COLLOCATED, qdata);
680 CeedOperatorSetField(op_apply, "v", Erestrictu, CEED_TRANSPOSE,
681 basisu, CEED_VECTOR_ACTIVE);
682
683 // Create the error operator
684 CeedOperatorCreate(ceed, qf_error, NULL, NULL, &op_error);
685 CeedOperatorSetField(op_error, "u", Erestrictu, CEED_TRANSPOSE,
686 basisu, CEED_VECTOR_ACTIVE);
687 CeedOperatorSetField(op_error, "true_soln", Erestrictui, CEED_NOTRANSPOSE,

--- 10 unchanged lines hidden (view full) ---

698 user->ltog0 = ltog0;
699 user->gtogD = gtogD;
700 }
701 user->Xloc = Xloc;
702 ierr = VecDuplicate(Xloc, &user->Yloc); CHKERRQ(ierr);
703 CeedVectorCreate(ceed, lsize*ncompu, &user->xceed);
704 CeedVectorCreate(ceed, lsize*ncompu, &user->yceed);
705 user->op = op_apply;
680 CeedOperatorSetField(op_apply, "v", Erestrictu, CEED_TRANSPOSE,
681 basisu, CEED_VECTOR_ACTIVE);
682
683 // Create the error operator
684 CeedOperatorCreate(ceed, qf_error, NULL, NULL, &op_error);
685 CeedOperatorSetField(op_error, "u", Erestrictu, CEED_TRANSPOSE,
686 basisu, CEED_VECTOR_ACTIVE);
687 CeedOperatorSetField(op_error, "true_soln", Erestrictui, CEED_NOTRANSPOSE,

--- 10 unchanged lines hidden (view full) ---

698 user->ltog0 = ltog0;
699 user->gtogD = gtogD;
700 }
701 user->Xloc = Xloc;
702 ierr = VecDuplicate(Xloc, &user->Yloc); CHKERRQ(ierr);
703 CeedVectorCreate(ceed, lsize*ncompu, &user->xceed);
704 CeedVectorCreate(ceed, lsize*ncompu, &user->yceed);
705 user->op = op_apply;
706 user->rho = rho;
706 user->qdata = qdata;
707 user->ceed = ceed;
708
709 ierr = MatCreateShell(comm, mnodes[0]*mnodes[1]*mnodes[2]*ncompu,
710 mnodes[0]*mnodes[1]*mnodes[2]*ncompu,
711 PETSC_DECIDE, PETSC_DECIDE, user, &mat); CHKERRQ(ierr);
712 if (bpChoice == CEED_BP1 || bpChoice == CEED_BP2) {
713 ierr = MatShellSetOperation(mat, MATOP_MULT, (void(*)(void))MatMult_Mass);
714 CHKERRQ(ierr);

--- 4 unchanged lines hidden (view full) ---

719 ierr = MatCreateVecs(mat, &rhs, NULL); CHKERRQ(ierr);
720
721 // Get RHS vector
722 ierr = VecDuplicate(Xloc, &rhsloc); CHKERRQ(ierr);
723 ierr = VecZeroEntries(rhsloc); CHKERRQ(ierr);
724 ierr = VecGetArray(rhsloc, &r); CHKERRQ(ierr);
725 CeedVectorSetArray(rhsceed, CEED_MEM_HOST, CEED_USE_POINTER, r);
726
707 user->ceed = ceed;
708
709 ierr = MatCreateShell(comm, mnodes[0]*mnodes[1]*mnodes[2]*ncompu,
710 mnodes[0]*mnodes[1]*mnodes[2]*ncompu,
711 PETSC_DECIDE, PETSC_DECIDE, user, &mat); CHKERRQ(ierr);
712 if (bpChoice == CEED_BP1 || bpChoice == CEED_BP2) {
713 ierr = MatShellSetOperation(mat, MATOP_MULT, (void(*)(void))MatMult_Mass);
714 CHKERRQ(ierr);

--- 4 unchanged lines hidden (view full) ---

719 ierr = MatCreateVecs(mat, &rhs, NULL); CHKERRQ(ierr);
720
721 // Get RHS vector
722 ierr = VecDuplicate(Xloc, &rhsloc); CHKERRQ(ierr);
723 ierr = VecZeroEntries(rhsloc); CHKERRQ(ierr);
724 ierr = VecGetArray(rhsloc, &r); CHKERRQ(ierr);
725 CeedVectorSetArray(rhsceed, CEED_MEM_HOST, CEED_USE_POINTER, r);
726
727 // Setup rho, rhs, and target
728 CeedOperatorApply(op_setup, xcoord, rho, CEED_REQUEST_IMMEDIATE);
727 // Setup qdata, rhs, and target
728 CeedOperatorApply(op_setup, xcoord, qdata, CEED_REQUEST_IMMEDIATE);
729 ierr = CeedVectorSyncArray(rhsceed, CEED_MEM_HOST); CHKERRQ(ierr);
730 CeedVectorDestroy(&xcoord);
731
732 // Gather RHS
733 ierr = VecRestoreArray(rhsloc, &r); CHKERRQ(ierr);
734 ierr = VecZeroEntries(rhs); CHKERRQ(ierr);
735 ierr = VecScatterBegin(ltog, rhsloc, rhs, ADD_VALUES, SCATTER_FORWARD);
736 CHKERRQ(ierr);

--- 100 unchanged lines hidden (view full) ---

837 ierr = VecScatterDestroy(&ltog); CHKERRQ(ierr);
838 ierr = VecScatterDestroy(&ltog0); CHKERRQ(ierr);
839 ierr = VecScatterDestroy(&gtogD); CHKERRQ(ierr);
840 ierr = MatDestroy(&mat); CHKERRQ(ierr);
841 ierr = KSPDestroy(&ksp); CHKERRQ(ierr);
842
843 CeedVectorDestroy(&user->xceed);
844 CeedVectorDestroy(&user->yceed);
729 ierr = CeedVectorSyncArray(rhsceed, CEED_MEM_HOST); CHKERRQ(ierr);
730 CeedVectorDestroy(&xcoord);
731
732 // Gather RHS
733 ierr = VecRestoreArray(rhsloc, &r); CHKERRQ(ierr);
734 ierr = VecZeroEntries(rhs); CHKERRQ(ierr);
735 ierr = VecScatterBegin(ltog, rhsloc, rhs, ADD_VALUES, SCATTER_FORWARD);
736 CHKERRQ(ierr);

--- 100 unchanged lines hidden (view full) ---

837 ierr = VecScatterDestroy(&ltog); CHKERRQ(ierr);
838 ierr = VecScatterDestroy(&ltog0); CHKERRQ(ierr);
839 ierr = VecScatterDestroy(&gtogD); CHKERRQ(ierr);
840 ierr = MatDestroy(&mat); CHKERRQ(ierr);
841 ierr = KSPDestroy(&ksp); CHKERRQ(ierr);
842
843 CeedVectorDestroy(&user->xceed);
844 CeedVectorDestroy(&user->yceed);
845 CeedVectorDestroy(&user->rho);
845 CeedVectorDestroy(&user->qdata);
846 CeedVectorDestroy(&target);
847 CeedOperatorDestroy(&op_setup);
848 CeedOperatorDestroy(&op_apply);
849 CeedOperatorDestroy(&op_error);
850 CeedElemRestrictionDestroy(&Erestrictu);
851 CeedElemRestrictionDestroy(&Erestrictx);
852 CeedElemRestrictionDestroy(&Erestrictui);
853 CeedElemRestrictionDestroy(&Erestrictxi);
854 CeedElemRestrictionDestroy(&Erestrictqdi);
855 CeedQFunctionDestroy(&qf_setup);
856 CeedQFunctionDestroy(&qf_apply);
857 CeedQFunctionDestroy(&qf_error);
858 CeedBasisDestroy(&basisu);
859 CeedBasisDestroy(&basisx);
860 CeedDestroy(&ceed);
861 ierr = PetscFree(user); CHKERRQ(ierr);
862 return PetscFinalize();
863}
846 CeedVectorDestroy(&target);
847 CeedOperatorDestroy(&op_setup);
848 CeedOperatorDestroy(&op_apply);
849 CeedOperatorDestroy(&op_error);
850 CeedElemRestrictionDestroy(&Erestrictu);
851 CeedElemRestrictionDestroy(&Erestrictx);
852 CeedElemRestrictionDestroy(&Erestrictui);
853 CeedElemRestrictionDestroy(&Erestrictxi);
854 CeedElemRestrictionDestroy(&Erestrictqdi);
855 CeedQFunctionDestroy(&qf_setup);
856 CeedQFunctionDestroy(&qf_apply);
857 CeedQFunctionDestroy(&qf_error);
858 CeedBasisDestroy(&basisu);
859 CeedBasisDestroy(&basisx);
860 CeedDestroy(&ceed);
861 ierr = PetscFree(user); CHKERRQ(ierr);
862 return PetscFinalize();
863}