[petsc-users] Problems with PetscFECreateLagrange(PETSC_COMM_WORLD, 3, 1, PETSC_TRUE, 2, 3, or 4, PETSC_DETERMINE, &fe_u)
Jed Brown
jed at jedbrown.org
Sat Jun 27 18:28:50 CDT 2026
- Previous message (by thread): [petsc-users] Problems with PetscFECreateLagrange(PETSC_COMM_WORLD, 3, 1, PETSC_TRUE, 2, 3, or 4, PETSC_DETERMINE, &fe_u)
- Next message (by thread): [petsc-users] Problems with PetscFECreateLagrange(PETSC_COMM_WORLD, 3, 1, PETSC_TRUE, 2, 3, or 4, PETSC_DETERMINE, &fe_u)
- Messages sorted by:
[ date ]
[ thread ]
[ subject ]
[ author ]
What do you mean by "with all the data obtained by the solver"?
One does not need any code to use this runtime option: -snes_view_solution vtk:solution.vtu
Are you wanting to store the sort of output you get from -snes_view into a file?
Rafel Amer Ramon via petsc-users <petsc-users at mcs.anl.gov> writes:
> Hi,
>
> I have a program ex06.c (I attach it) that uses PETSc and
>
> PetscFECreateLagrange(PETSC_COMM_WORLD, 3, 1, PETSC_TRUE, 2, 3, or 4, PETSC_DETERMINE, &fe_u)
>
> When I save the solution to a file solucio_torus.vtu with the code
>
> PetscViewer viewer_vtu;
> PetscSection section;
>
> ierr = PetscViewerCreate(PETSC_COMM_WORLD, &viewer_vtu); CHKERRQ(ierr);
> ierr = PetscViewerSetType(viewer_vtu, PETSCVIEWERVTK); CHKERRQ(ierr);
>
> ierr = PetscViewerFileSetMode(viewer_vtu, FILE_MODE_WRITE); CHKERRQ(ierr);
> ierr = PetscViewerFileSetName(viewer_vtu, "solucio_torus.vtu"); CHKERRQ(ierr);
>
> ierr = DMGetLocalSection(dm, §ion); CHKERRQ(ierr);
> if (section)
> ierr = PetscSectionSetFieldName(section, 0, "Solucio"); CHKERRQ(ierr);
>
> ierr = PetscObjectSetName((PetscObject)u_solucio, ""); CHKERRQ(ierr);
> ierr = VecView(u_solucio, viewer_vtu); CHKERRQ(ierr);
> ierr = PetscViewerDestroy(&viewer_vtu); CHKERRQ(ierr);
>
> it always generate a file solucio_torus.vtu of 4.5MB independently of the degree of the Lagrange polynomials, 2, 3 or 4.
> I run the code
>
> mpirun -np 32 --hostfile ~/hosts ./ex08 -mesh_input torus.msh -lua_input function2.lua -order 3 -snes_type ksponly -ksp_type gmres -pc_type gamg
>
> Does anybody knows how to save a file solucio_torus.vtu
> with all the data obtained by the solver
> ierr = SNESSolve(snes, NULL, u_solucio); CHKERRQ(ierr);
> Thank you very much.
>
> Best regards,
>
> Rafel Amer
> /**************************************************************************************
> * Filename: ex04.c
> * Author: Rafel Amer (rafel.amer AT upc.edu)
> * Copyright: Rafel Amer 2026
> * Disclaimer: This code is presented "as is" and it has been written for
> * educational purposes.
> *
> * License: This library is free software; you can redistribute it and/or
> * modify it under the terms of either:
> *
> * 1 the GNU Lesser General Public License as published by the Free
> * Software Foundation; either version 3 of the License, or (at your
> * option) any later version.
> *
> * or
> *
> * 2 the GNU General Public License as published by the Free Software
> * Foundation; either version 2 of the License, or (at your option)
> * any later version.
> *
> * See https://urldefense.us/v3/__https://www.gnu.org/licenses/__;!!G_uCfscf7eWS!e8Q8RJGC0Wq7jBFniRT82ho-dQ4iQN0Jzo4li6tbvUjy9c5XQZxVLmJ_XZiAx2pOKG7DrS9Q7k0-Qabt7Oo$
> ***************************************************************************************/
> #include <petscdmplex.h>
> #include <petscsys.h>
> #include <petscds.h>
> #include <petscfe.h>
> #include <petscsnes.h> // Canviem a la llibreria SNES per a la gestió correcta de memòria
> #include <petscviewerhdf5.h>
> #include <lua.h>
> #include <lualib.h>
> #include <lauxlib.h>
>
> static char help[] = "Solucionador del tòrus via SNES (ksponly) per a elements finits en paral·lel.\n\n";
>
> static PetscErrorCode CondicioContornDirichlet(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)
> {
> // u[0] = 0.0;
> u[0] = 100*x[2];
> return 0;
> }
>
> static void f0_volum(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[])
> {
> PetscScalar c = (a && NfAux > 0 && aOff) ? a[aOff[0] + 4] : 0.0;
> PetscScalar f = (a && NfAux > 0 && aOff) ? a[aOff[0] + 5] : 0.0;
> f0[0] = c * u[0] - f;
> }
>
> static void f1_volum(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[])
> {
> PetscScalar coef_a = (a && NfAux > 0 && aOff) ? a[aOff[0] + 0] : 1.0;
> PetscScalar b1 = (a && NfAux > 0 && aOff) ? a[aOff[0] + 1] : 0.0;
> PetscScalar b2 = (a && NfAux > 0 && aOff) ? a[aOff[0] + 2] : 0.0;
> PetscScalar b3 = (a && NfAux > 0 && aOff) ? a[aOff[0] + 3] : 0.0;
>
> f1[0] = coef_a * u_x[0] + b1 * u[0];
> f1[1] = coef_a * u_x[1] + b2 * u[0];
> f1[2] = coef_a * u_x[2] + b3 * u[0];
> }
>
> static void g3_jacobia(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, PetscReal u_tShift, const PetscReal x[], PetscInt numConstants, const PetscScalar constants[], PetscScalar g3[])
> {
> PetscScalar coef_a = (a && NfAux > 0 && aOff) ? a[aOff[0] + 0] : 1.0;
> g3[0] = coef_a; g3[1] = 0.0; g3[2] = 0.0;
> g3[3] = 0.0; g3[4] = coef_a; g3[5] = 0.0;
> g3[6] = 0.0; g3[7] = 0.0; g3[8] = coef_a;
> }
>
> static PetscErrorCode EvaluaCoeficientsLua(PetscInt dim, PetscReal time, const PetscReal x[], PetscInt Nc, PetscScalar u[], void *ctx)
> {
> lua_State *L = (lua_State *)ctx;
> lua_getglobal(L, "fem");
> lua_pushnumber(L, x[0]);
> lua_pushnumber(L, x[1]);
> lua_pushnumber(L, x[2]);
> if (lua_pcall(L, 3, 1, 0) != LUA_OK) return PETSC_ERR_USER;
>
> lua_getfield(L, -1, "a"); u[0] = lua_tonumber(L, -1); lua_pop(L, 1);
> lua_getfield(L, -1, "b1"); u[1] = lua_tonumber(L, -1); lua_pop(L, 1);
> lua_getfield(L, -1, "b2"); u[2] = lua_tonumber(L, -1); lua_pop(L, 1);
> lua_getfield(L, -1, "b3"); u[3] = lua_tonumber(L, -1); lua_pop(L, 1);
> lua_getfield(L, -1, "c"); u[4] = lua_tonumber(L, -1); lua_pop(L, 1);
> lua_getfield(L, -1, "f"); u[5] = lua_tonumber(L, -1); lua_pop(L, 1);
> lua_pop(L, 1);
> return 0;
> }
>
> int main(int argc, char **argv)
> {
> DM dm, dmDist = NULL;
> PetscErrorCode ierr;
> PetscMPIInt rank, size;
> PetscBool flg;
> lua_State *L;
> PetscInt ordre_element = 4;
> char filename[PETSC_MAX_PATH_LEN] = "torus.msh";
> char luafilename[PETSC_MAX_PATH_LEN] = "function.lua";
>
> ierr = PetscInitialize(&argc, &argv, (char*)0, help); if (ierr) return ierr;
> ierr = MPI_Comm_rank(PETSC_COMM_WORLD, &rank); CHKERRQ(ierr);
> ierr = MPI_Comm_size(PETSC_COMM_WORLD, &size); CHKERRQ(ierr);
>
> ierr = PetscOptionsGetString(NULL, NULL, "-mesh_input", filename, sizeof(filename), &flg); CHKERRQ(ierr);
> ierr = PetscOptionsGetString(NULL, NULL, "-lua_input", luafilename, sizeof(luafilename), &flg); CHKERRQ(ierr);
> ierr = PetscOptionsGetInt(NULL, NULL, "-order", &ordre_element, &flg); CHKERRQ(ierr);
>
> // [Pas 1] Llegir i distribuir la malla
> ierr = DMPlexCreateGmshFromFile(PETSC_COMM_WORLD, filename, PETSC_TRUE, &dm); CHKERRQ(ierr);
> ierr = DMSetFromOptions(dm); CHKERRQ(ierr);
>
> // Forcem a PETSc a incloure les etiquetes de Gmsh durant la distribució paral·lela
> ierr = DMPlexDistribute(dm, 0, NULL, &dmDist);
> CHKERRQ(ierr);
> if (dmDist)
> {
> ierr = DMDestroy(&dm); CHKERRQ(ierr);
> dm = dmDist;
> }
>
> // [Pas 2] Carregar Lua de forma aïllada
> if ((L = luaL_newstate()) == NULL) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_MEM, "Error creant Lua.");
> luaL_openlibs(L);
> if (luaL_dofile(L, luafilename) != LUA_OK) SETERRQ(PETSC_COMM_SELF, PETSC_ERR_FILE_OPEN, "Error carregant function.lua");
>
> // [Pas 3] Configurar l'espai discret principal (u)
> PetscFE fe_u;
> ierr = PetscFECreateLagrange(PETSC_COMM_WORLD, 3, 1, PETSC_TRUE, ordre_element, PETSC_DETERMINE, &fe_u); CHKERRQ(ierr);
> ierr = DMSetField(dm, 0, NULL, (PetscObject)fe_u); CHKERRQ(ierr);
> ierr = DMCreateDS(dm); CHKERRQ(ierr);
>
> PetscDS ds;
> ierr = DMGetDS(dm, &ds); CHKERRQ(ierr);
> ierr = PetscDSSetResidual(ds, 0, f0_volum, f1_volum); CHKERRQ(ierr);
> ierr = PetscDSSetJacobian(ds, 0, 0, NULL, NULL, NULL, g3_jacobia); CHKERRQ(ierr);
>
> DMLabel label;
> // 1. Intentem buscar el nom que has posat tu
> ierr = DMGetLabel(dm, "Frontera", &label); CHKERRQ(ierr);
> if (label && rank == 0) PetscPrintf(PETSC_COMM_WORLD, "--> Trobat per nom: 'Frontera'\n");
>
> // 2. Si no, usem el marker estàndard de PETSc
> if (!label) {
> ierr = DMGetLabel(dm, "marker", &label); CHKERRQ(ierr);
> if (label && rank == 0) PetscPrintf(PETSC_COMM_WORLD, "--> Trobat per nom genèric: 'marker'\n");
> }
> if (!label) { ierr = DMGetLabel(dm, "Face Sets", &label); CHKERRQ(ierr); }
>
> // MODIFICACIÓ DE SEGURETAT: Apliquem a l'ID 1 i a l'ID 2 simultàniament
> // Si Gmsh o PETSc han mogut l'ID de la superfície del 1 al 2 durant la conversió, els cobrim tots dos.
> PetscInt numIds = 2;
> PetscInt ids[2] = {1, 2};
>
> if (label) {
> ierr = PetscDSAddBoundary(ds, DM_BC_ESSENTIAL, "ContornExterior", label,
> numIds, ids, 0, 0, NULL,
> (PetscVoidFn *)CondicioContornDirichlet, NULL, NULL, NULL);
> CHKERRQ(ierr);
> } else {
> PetscPrintf(PETSC_COMM_WORLD, "ALERTA: No s'ha trobat cap etiqueta de contorn!\n");
> }
>
> // ====================================================================
> // [Pas 4] ESTRATÈGIA EX03 SENSE INTERFERÈNCIES DE LUA
> // ====================================================================
> DM dmAux;
> PetscFE fe_aux;
> Vec vectorTemporalLua, vectorResultatsFinal;
>
> ierr = DMClone(dm, &dmAux); CHKERRQ(ierr);
> ierr = PetscFECreateLagrange(PETSC_COMM_WORLD, 3, 6, PETSC_TRUE, ordre_element, PETSC_DETERMINE, &fe_aux); CHKERRQ(ierr);
> ierr = DMSetField(dmAux, 0, NULL, (PetscObject)fe_aux); CHKERRQ(ierr);
> ierr = DMCreateDS(dmAux); CHKERRQ(ierr);
>
> ierr = DMCreateLocalVector(dmAux, &vectorTemporalLua); CHKERRQ(ierr);
>
> PetscErrorCode (*funcs[1])(PetscInt, PetscReal, const PetscReal[], PetscInt, PetscScalar[], void *);
> void *ctxs[1];
> funcs[0] = EvaluaCoeficientsLua;
> ctxs[0] = L;
>
> ierr = DMProjectFunctionLocal(dmAux, 0.0, funcs, ctxs, INSERT_VALUES, vectorTemporalLua); CHKERRQ(ierr);
>
> // Tanquem LUA immediatament, exactament com feies a l'ex03 per no embrutar la memòria
> lua_gc(L, LUA_GCCOLLECT, 0);
> lua_close(L);
>
> // Dupliquem per assegurar-nos un vector completament net i homogeni de PETSc
> ierr = DMCreateLocalVector(dmAux, &vectorResultatsFinal); CHKERRQ(ierr);
> ierr = VecCopy(vectorTemporalLua, vectorResultatsFinal); CHKERRQ(ierr);
> ierr = VecDestroy(&vectorTemporalLua); CHKERRQ(ierr);
>
> ierr = DMSetAuxiliaryVec(dm, NULL, 0, 0, vectorResultatsFinal); CHKERRQ(ierr);
>
> // [Pas 5] Creem el vector solució global
> Vec u_solucio;
> ierr = DMCreateGlobalVector(dm, &u_solucio); CHKERRQ(ierr);
> ierr = VecZeroEntries(u_solucio); CHKERRQ(ierr);
>
> // ====================================================================
> // [Pas 6] CONFIGURACIÓ DE SNES (EL MOTOR CORRECTE PER A DMPLEX FEM)
> // ====================================================================
> SNES snes;
> ierr = SNESCreate(PETSC_COMM_WORLD, &snes); CHKERRQ(ierr);
> ierr = SNESSetDM(snes, dm); CHKERRQ(ierr);
>
> // Establir les funcions de DMPlex per al residu i jacobià reals
> ierr = DMSNESSetFunctionLocal(dm, (PetscErrorCode (*)(DM, Vec, Vec, void *))DMPlexSNESComputeResidualFEM, NULL); CHKERRQ(ierr);
> ierr = DMSNESSetJacobianLocal(dm, (PetscErrorCode (*)(DM, Vec, Mat, Mat, void *))DMPlexSNESComputeJacobianFEM, NULL); CHKERRQ(ierr);
>
> ierr = SNESSetFromOptions(snes); CHKERRQ(ierr);
> // Forcem el solver a comportar-se com un KSP lineal pur
> ierr = SNESSetType(snes, SNESKSPONLY); CHKERRQ(ierr);
>
> //-----------------------------------------------------------------
> // ORDENAT CORRECTAMENT: Declarar, extreure i llavors configurar
> //-----------------------------------------------------------------
> KSP ksp; // 1. Declarem la variable
> ierr = SNESGetKSP(snes, &ksp); CHKERRQ(ierr); // 2. L'extraiem del SNES
> ierr = KSPSetType(ksp, KSPGMRES); CHKERRQ(ierr); // 3. Configurem el tipus (GMRES)
> //-----------------------------------------------------------------
>
> // Executem la resolució paral·lela
> ierr = SNESSolve(snes, NULL, u_solucio); CHKERRQ(ierr);
>
> // ===========================================================
> // PAS 7: ESCRIPTURA NATIVA EN FORMAT VTK (.vtu) - NET DEFINITIU
> // ===========================================================
> PetscViewer viewer_vtu;
> PetscSection section;
>
> // 1. Creem el viewer en format VTK pur
> ierr = PetscViewerCreate(PETSC_COMM_WORLD, &viewer_vtu); CHKERRQ(ierr);
> ierr = PetscViewerSetType(viewer_vtu, PETSCVIEWERVTK); CHKERRQ(ierr);
>
> // 2. Configurem el fitxer de sortida
> ierr = PetscViewerFileSetMode(viewer_vtu, FILE_MODE_WRITE); CHKERRQ(ierr);
> ierr = PetscViewerFileSetName(viewer_vtu, "solucio_torus.vtu"); CHKERRQ(ierr);
>
> // 3. LA MÀGIA: Definim el nom del camp net a la secció
> ierr = DMGetLocalSection(dm, §ion); CHKERRQ(ierr);
> if (section) {
> ierr = PetscSectionSetFieldName(section, 0, "Solucio"); CHKERRQ(ierr);
> }
>
> // 4. EL TRUC: No li donem el mateix nom al vector (deixem-lo buit o genèric)
> // d'aquesta manera evitem la duplicació "SolucioSolucio"
> ierr = PetscObjectSetName((PetscObject)u_solucio, ""); CHKERRQ(ierr);
>
> // 5. Escriptura del vector interpolat pels nodes de la malla
> ierr = VecView(u_solucio, viewer_vtu); CHKERRQ(ierr);
>
> // 6. Destruïm el viewer
> ierr = PetscViewerDestroy(&viewer_vtu); CHKERRQ(ierr);
>
> // ====================================================================
> // [Pas 8] ESTADÍSTIQUES FINALS
> // ====================================================================
> PetscInt its;
> SNESConvergedReason reason;
> ierr = SNESGetIterationNumber(snes, &its); CHKERRQ(ierr);
> ierr = SNESGetConvergedReason(snes, &reason); CHKERRQ(ierr);
>
> PetscInt n_equacions;
> ierr = VecGetSize(u_solucio, &n_equacions); CHKERRQ(ierr);
>
> if (rank == 0) {
> PetscPrintf(PETSC_COMM_WORLD, "===========================================================\n");
> PetscPrintf(PETSC_COMM_WORLD, " RESOLUCIÓ PARAL·LELA DE L'EQUACIÓ DEL TÒRUS AMB ÈXIT\n");
> PetscPrintf(PETSC_COMM_WORLD, " Graus de llibertat totals a la malla: %d\n", n_equacions);
> PetscPrintf(PETSC_COMM_WORLD, " Raó de convergència del Solver (Reason ID): %d\n", reason);
> PetscPrintf(PETSC_COMM_WORLD, " Iteracions realitzades: %d\n", its);
> PetscPrintf(PETSC_COMM_WORLD, "===========================================================\n");
>
> PetscPrintf(PETSC_COMM_WORLD, "===========================================================\n");
> PetscPrintf(PETSC_COMM_WORLD, " Generat el fitxer pur 'solucio_torus.vtu'\n");
> PetscPrintf(PETSC_COMM_WORLD, "===========================================================\n");
> }
> // Alliberament final impecable
> ierr = SNESDestroy(&snes); CHKERRQ(ierr);
> ierr = VecDestroy(&u_solucio); CHKERRQ(ierr);
> ierr = VecDestroy(&vectorResultatsFinal); CHKERRQ(ierr);
> ierr = DMDestroy(&dmAux); CHKERRQ(ierr);
> ierr = PetscFEDestroy(&fe_u); CHKERRQ(ierr);
> ierr = PetscFEDestroy(&fe_aux); CHKERRQ(ierr);
> ierr = DMDestroy(&dm);
> ierr = PetscFinalize();
> return ierr;
> }
> function a(x,y,z)
> return 1.0 -- Sempre serà >= 1.0, mai zero!
> end
>
> function b1(x,y,z) return 0.0 end -- Un transport molt suau
> function b2(x,y,z) return 0.0 end
> function b3(x,y,z) return 0.0 end
>
> function c(x,y,z)
> return 0.0 -- Reacció amortidora
> end
>
> function f(x,y,z)
> return 1.0*z
> end
>
> function fem(x,y,z)
> local r = {}
> r["a"] = a(x,y,z)
> r["b1"] = b1(x,y,z)
> r["b2"] = b2(x,y,z)
> r["b3"] = b3(x,y,z)
> r["c"] = c(x,y,z)
> r["f"] = f(x,y,z)
> return r
> end
>
> function funcio_contorn(x,y,z)
> return 100*z
> end
- Previous message (by thread): [petsc-users] Problems with PetscFECreateLagrange(PETSC_COMM_WORLD, 3, 1, PETSC_TRUE, 2, 3, or 4, PETSC_DETERMINE, &fe_u)
- Next message (by thread): [petsc-users] Problems with PetscFECreateLagrange(PETSC_COMM_WORLD, 3, 1, PETSC_TRUE, 2, 3, or 4, PETSC_DETERMINE, &fe_u)
- Messages sorted by:
[ date ]
[ thread ]
[ subject ]
[ author ]
More information about the petsc-users
mailing list