v0.10.0
error_norm.cpp
Go to the documentation of this file.
1 /*
2  \file error_norm.cpp
3  \ingroup mofem_helmholtz_elem
4 
5  Calculates finite element (Galerkin) approximation for difference between two solutions in L^2 and H_1 norm.
6  \bug work for scalar filed only, NEED further modification.
7  */
8 
9 /*
10  * This file is part of MoFEM.
11  * MoFEM is free software: you can redistribute it and/or modify it under
12  * the terms of the GNU Lesser General Public License as published by the
13  * Free Software Foundation, either version 3 of the License, or (at your
14  * option) any later version.
15  *
16  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
17  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
18  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
19  * License for more details.
20  *
21  * You should have received a copy of the GNU Lesser General Public
22  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>.
23  *
24  */
25 
26 #include <MoFEM.hpp>
27 using namespace MoFEM;
28 #include <PostProcOnRefMesh.hpp>
29 #include <NormElement.hpp>
30 
32 #include <boost/numeric/ublas/vector_proxy.hpp>
33 #include <petsctime.h>
34 #include <fstream>
35 #include <iostream>
36 #include <stdexcept>
37 
38 static char help[] = "...\n\n";
39 
40 int main(int argc, char *argv[]) {
41 
42  ErrorCode rval;
43  PetscErrorCode ierr;
44 
45  PetscInitialize(&argc,&argv,(char *)0,help);
46 
47  moab::Core mb_instance;
48  moab::Interface& moab = mb_instance;
49  int rank;
50  MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
51 
52  //read h5m solution file into moab
53  bool usel2; //norm type
54  PetscBool flg = PETSC_TRUE;
55 
56  char mesh_file_name[255];
57  ierr = PetscOptionsGetString(PETSC_NULL,"-my_file",mesh_file_name,255,&flg); CHKERRQ(ierr);
58  if(flg != PETSC_TRUE) {
59  SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -my_file (MESH FILE NEEDED)");
60  }
61 
62 
63  //PetscBool userela = PETSC_FALSE;
64  //ierr = PetscOptionsGetBool(NULL,"-relative_error",&userela,NULL); CHKERRQ(ierr);
65 
66  /* FIX ME, change to enum for better performance */
67  char type_error_norm[255];
68  ierr = PetscOptionsGetString(PETSC_NULL,"-norm_type",type_error_norm,255,&flg); CHKERRQ(ierr);
69  if(flg != PETSC_TRUE) {
70  SETERRQ(PETSC_COMM_SELF,1,"*** ERROR -type_error_norm (L2 or H1 type needed)");
71  }
72 
73  if (strcmp ("l2",type_error_norm ) == 0) {usel2 = true;}
74  else if(strcmp ("h1",type_error_norm ) == 0) {usel2 = false;}
75 
76  ParallelComm* pcomm = ParallelComm::get_pcomm(&moab,MYPCOMM_INDEX);
77  if(pcomm == NULL) pcomm = new ParallelComm(&moab,PETSC_COMM_WORLD);
78 
79  const char *option;
80  option = "";//"PARALLEL=BCAST;";//;DEBUG_IO";
81  // BARRIER_RANK_START(pcomm)
82  /* load the mesh files */
83  rval = moab.load_file(mesh_file_name, 0, option); CHKERRQ_MOAB(rval);
84  //rval = moab.load_file(mesh_file_name2, 0, option); CHKERRQ_MOAB(rval);
85  // BARRIER_RANK_END(pcomm)
86 
87  //Create MoFEM (Joseph) database
88  MoFEM::Core core(moab);
89  MoFEM::Interface& m_field = core;
90 
91  //count the comsumption of time by single run
92  PetscLogDouble t1,t2;
93  PetscLogDouble v1,v2;
94  ierr = PetscTime(&v1); CHKERRQ(ierr);
95  ierr = PetscGetCPUTime(&t1); CHKERRQ(ierr);
96 
97  //set entitities bit level
98  BitRefLevel bit_level0;
99  bit_level0.set(0);
100  EntityHandle meshset_level0;
101  rval = moab.create_meshset(MESHSET_SET,meshset_level0); CHKERRQ_MOAB(rval);
102  ierr = m_field.seed_ref_level_3D(0,bit_level0); CHKERRQ(ierr);
103 
104  //Fields
105  ierr = m_field.add_field("erorNORM",H1,AINSWORTH_LOBATTO_BASE,1,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
106 
107  //Problem
108  ierr = m_field.add_problem("NORM_PROBLEM"); CHKERRQ(ierr);
109 
110  //set refinment level for problem
111  ierr = m_field.modify_problem_ref_level_add_bit("NORM_PROBLEM",bit_level0); CHKERRQ(ierr);
112 
113 
114 
115 
116  //meshset consisting all entities in mesh
117  EntityHandle root_set = moab.get_root_set();
118  //add entities to field
119  ierr = m_field.add_ents_to_field_by_TETs(root_set,"erorNORM"); CHKERRQ(ierr);
120 
121 
122  //set app. order , approximation error of norm should be as least 1 order higher than numerical spaces.
123  //see Hierarchic Finite Element Bases on Unstructured Tetrahedral Meshes (Mark Ainsworth & Joe Coyle)
124  PetscInt order;
125  ierr = PetscOptionsGetInt(PETSC_NULL,"-my_order",&order,&flg); CHKERRQ(ierr);
126  if(flg != PETSC_TRUE) {
127  order = 2;
128  }
129 
130  ierr = m_field.set_field_order(root_set,MBTET,"erorNORM",order); CHKERRQ(ierr);
131  ierr = m_field.set_field_order(root_set,MBTRI,"erorNORM",order); CHKERRQ(ierr);
132  ierr = m_field.set_field_order(root_set,MBEDGE,"erorNORM",order); CHKERRQ(ierr);
133  ierr = m_field.set_field_order(root_set,MBVERTEX,"erorNORM",1); CHKERRQ(ierr);
134 
135  if(!m_field.check_field("MESH_NODE_POSITIONS")) {
136  ierr = m_field.add_field("MESH_NODE_POSITIONS",H1,AINSWORTH_LEGENDRE_BASE,3,MB_TAG_SPARSE,MF_ZERO); CHKERRQ(ierr);
137  ierr = m_field.add_ents_to_field_by_TETs(root_set,"MESH_NODE_POSITIONS"); CHKERRQ(ierr);
138  ierr = m_field.set_field_order(0,MBTET,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
139  ierr = m_field.set_field_order(0,MBTRI,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
140  ierr = m_field.set_field_order(0,MBEDGE,"MESH_NODE_POSITIONS",2); CHKERRQ(ierr);
141  ierr = m_field.set_field_order(0,MBVERTEX,"MESH_NODE_POSITIONS",1); CHKERRQ(ierr);
142  }
143 
144  /* global error and analytical solution */
145  double error;
146  double analy;
147 
148 
149  NormElement norm_elements(m_field,error,analy);
150 
151  if(m_field.check_field("reEX") && m_field.check_field("rePRES")) {
152  //&& m_field.check_field("imPRES") && m_field.check_field("imEX") ) {
153  norm_elements.addNormElements("NORM_PROBLEM","NORM_FE","erorNORM","reEX","rePRES","imEX","imPRES");
154  }
155 
156 
157  /****/
158  //build database
159 
160  //build field
161  ierr = m_field.build_fields(); CHKERRQ(ierr);
162  //build finite elemnts
163  ierr = m_field.build_finite_elements(); CHKERRQ(ierr);
164  //build adjacencies
165  ierr = m_field.build_adjacencies(bit_level0); CHKERRQ(ierr);
166  //build problem
167  // ierr = m_field.build_problems(); CHKERRQ(ierr);
168  ierr = m_field.build_problem("NORM_PROBLEM",true); CHKERRQ(ierr);
169 
170  Projection10NodeCoordsOnField ent_method_material(m_field,"MESH_NODE_POSITIONS");
171  ierr = m_field.loop_dofs("MESH_NODE_POSITIONS",ent_method_material); CHKERRQ(ierr);
172 
173 
174  /****/
175  //mesh partitioning
176  //partition
177  ierr = m_field.partition_problem("NORM_PROBLEM"); CHKERRQ(ierr);
178  ierr = m_field.partition_finite_elements("NORM_PROBLEM"); CHKERRQ(ierr);
179  //what are ghost nodes, see Petsc Manual
180  ierr = m_field.partition_ghost_dofs("NORM_PROBLEM"); CHKERRQ(ierr);
181 
182  if(m_field.check_field("reEX") && m_field.check_field("imEX")) {
183 
184  ierr = m_field.build_problem("EX1_PROBLEM",true); CHKERRQ(ierr);
185  ierr = m_field.partition_problem("EX1_PROBLEM"); CHKERRQ(ierr);
186  ierr = m_field.partition_finite_elements("EX1_PROBLEM"); CHKERRQ(ierr);
187  ierr = m_field.partition_ghost_dofs("EX1_PROBLEM"); CHKERRQ(ierr);
188 
189  }
190 
191  //print block sets with materials
192  //ierr = m_field.print_cubit_materials_set(); CHKERRQ(ierr);
193 
194  Vec F;
195  ierr = m_field.VecCreateGhost("NORM_PROBLEM",ROW,&F); CHKERRQ(ierr);
196  Vec T;
197  ierr = VecDuplicate(F,&T); CHKERRQ(ierr);
198  Mat A;
199  ierr = m_field.MatCreateMPIAIJWithArrays("NORM_PROBLEM",&A); CHKERRQ(ierr);
200 
201 
202  /* Set operators */
203  ierr = norm_elements.setNormFiniteElementRhsOperator("erorNORM","reEX","rePRES","imEX","imPRES",F,usel2); CHKERRQ(ierr);
204  ierr = norm_elements.setNormFiniteElementLhsOperator("erorNORM",A); CHKERRQ(ierr);
205 
206 
207 
208  /* create matrix and vectors */
209  ierr = VecZeroEntries(T); CHKERRQ(ierr);
210  ierr = VecGhostUpdateBegin(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
211  ierr = VecGhostUpdateEnd(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
212  ierr = VecZeroEntries(F); CHKERRQ(ierr);
213  ierr = VecGhostUpdateBegin(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
214  ierr = VecGhostUpdateEnd(F,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
215  ierr = MatZeroEntries(A); CHKERRQ(ierr);
216  ierr = m_field.set_global_ghost_vector("NORM_PROBLEM",ROW,T,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
217 
218  /* looop over finite elements */
219  ierr = m_field.loop_finite_elements("NORM_PROBLEM","NORM_FE",norm_elements.getLoopFe()); CHKERRQ(ierr);
220 
221 
222  /* create ghost points */
223  ierr = VecGhostUpdateBegin(F,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
224  ierr = VecGhostUpdateEnd(F,ADD_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
225  ierr = VecAssemblyBegin(F); CHKERRQ(ierr);
226  ierr = VecAssemblyEnd(F); CHKERRQ(ierr);
227  ierr = MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
228  ierr = MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY); CHKERRQ(ierr);
229 
230  //Solver
231  KSP solver1;
232  ierr = KSPCreate(PETSC_COMM_WORLD,&solver1); CHKERRQ(ierr);
233  ierr = KSPSetOperators(solver1,A,A); CHKERRQ(ierr);
234  ierr = KSPSetFromOptions(solver1); CHKERRQ(ierr);
235  ierr = KSPSetUp(solver1); CHKERRQ(ierr);
236 
237  ierr = KSPSolve(solver1,F,T); CHKERRQ(ierr);
238  ierr = VecGhostUpdateBegin(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
239  ierr = VecGhostUpdateEnd(T,INSERT_VALUES,SCATTER_FORWARD); CHKERRQ(ierr);
240  ierr = m_field.set_global_ghost_vector("NORM_PROBLEM",ROW,T,INSERT_VALUES,SCATTER_REVERSE); CHKERRQ(ierr);
241 
242  /* Global error calculation */
243 
244  //PetscPrintf(PETSC_COMM_WORLD,"\n real part of l2 norm of analytical solution is: %f \n",real_analy);
245  //PetscPrintf(PETSC_COMM_WORLD,"\n imag part of l2 norm of analytical solution is: %f \n",imag_analy);
246 
247  if(usel2) {
248  double aa = sqrt(error)/sqrt(analy);
249  PetscPrintf(PETSC_COMM_WORLD,"\n global l2 realtive error is: %f \n",aa);
250  if(!pcomm->rank()) {
251  std::cout << "precise l2 realtive error = "<< std::setprecision(16) << aa << '\n'; /* 13 digits precision */
252  }
253  } else {
254 
255  double aa = sqrt(error)/sqrt(analy);
256  PetscPrintf(PETSC_COMM_WORLD,"\n global H1 realtive error is: %f \n",aa);
257  if(!pcomm->rank()) {
258  std::cout << "precise h1 realtive error = "<< std::setprecision(16) << aa << '\n'; /* 13 digits precision */
259  }
260  }
261 
262  /* destroy objects */
263  ierr = MatDestroy(&A); CHKERRQ(ierr);
264  ierr = VecDestroy(&F); CHKERRQ(ierr);
265  ierr = VecDestroy(&T); CHKERRQ(ierr);
266  ierr = KSPDestroy(&solver1); CHKERRQ(ierr);
267 
268  /* save local error indicator on mesh */
269  PetscBool save_postproc_mesh = PETSC_TRUE;
270  ierr = PetscOptionsGetBool(NULL,"-save_postproc_mesh",&save_postproc_mesh,NULL); CHKERRQ(ierr);
271  if(save_postproc_mesh) {
272 
273  PostProcVolumeOnRefinedMesh post_proc1(m_field);
275  ierr = post_proc1.addFieldValuesPostProc("erorNORM"); CHKERRQ(ierr);
276  ierr = post_proc1.addFieldValuesPostProc("MESH_NODE_POSITIONS"); CHKERRQ(ierr);
277  ierr = m_field.loop_finite_elements("NORM_PROBLEM","NORM_FE",post_proc1); CHKERRQ(ierr);
278  ierr = post_proc1.writeFile("norm_error.h5m"); CHKERRQ(ierr);
279 
280  }
281 
282  ierr = PetscTime(&v2);CHKERRQ(ierr);
283  ierr = PetscGetCPUTime(&t2);CHKERRQ(ierr);
284  PetscSynchronizedPrintf(PETSC_COMM_WORLD,"Total Rank %d Time = %f S CPU Time = %f S \n",pcomm->rank(),v2-v1,t2-t1);
285 
286  ierr = PetscFinalize(); CHKERRQ(ierr);
287 
288  return 0;
289 
290 }
#define CHKERRQ_MOAB(a)
check error code of MoAB function
Definition: definitions.h:523
const double T
finite element to appeximate analytical solution on surface
Deprecated interface functions.
virtual MoFEMErrorCode build_adjacencies(const Range &ents, int verb=DEFAULT_VERBOSITY)=0
build adjacencies
DEPRECATED MoFEMErrorCode set_global_ghost_vector(const Problem *problem_ptr, RowColData rc, Vec V, InsertMode mode, ScatterMode scatter_mode) const
set values of vector from/to mesh database (collective)collective - need tu be run on all processors ...
CoreTmp< 0 > Core
Definition: Core.hpp:1129
Post-process fields on refined mesh.
static char help[]
Definition: error_norm.cpp:38
virtual MoFEMErrorCode add_field(const std::string &name, const FieldSpace space, const FieldApproximationBase base, const FieldCoefficientsNumber nb_of_coefficients, const TagType tag_type=MB_TAG_SPARSE, const enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add field.
virtual MoFEMErrorCode set_field_order(const EntityHandle meshset, const EntityType type, const std::string &name, const ApproximationOrder order, int verb=DEFAULT_VERBOSITY)=0
Set order approximation of the entities in the field.
Operators and data structures for L^2Norm analysis.
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
char mesh_file_name[255]
Projection of edge entities with one mid-node on hierarchical basis.
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:85
MoFEMErrorCode writeFile(const std::string file_name)
wrote results in (MOAB) format, use "file_name.h5m"
PetscErrorCode PetscOptionsGetString(PetscOptions *, const char pre[], const char name[], char str[], size_t size, PetscBool *set)
DEPRECATED MoFEMErrorCode partition_finite_elements(const std::string &name, bool part_from_moab=false, int low_proc=-1, int hi_proc=-1, int verb=-1)
partition finite elements
MoFEMErrorCode addFieldValuesPostProc(const std::string field_name, Vec v=PETSC_NULL)
Add operator to post-process L2, H1, Hdiv, Hcurl field value.
PetscErrorCode addNormElements(const string problem, string fe, const string norm_field_name, const string field1_name, const string field2_name, const string mesh_nodals_positions="MESH_NODE_POSITIONS")
constexpr int order
CHKERRQ(ierr)
int main(int argc, char *argv[])
Definition: error_norm.cpp:40
Ainsworth Cole (Legendre) approx. base .
Definition: definitions.h:152
DEPRECATED MoFEMErrorCode build_problem(const std::string &name, const bool square_matrix, int verb=-1)
build problem data structures
PetscErrorCode setNormFiniteElementRhsOperator(string norm_field_name, string field1_name, string field2_name, Mat A, Vec &F, bool usel2, bool userela, string nodals_positions="MESH_NODE_POSITIONS")
Post processing.
virtual MoFEMErrorCode loop_dofs(const Problem *problem_ptr, const std::string &field_name, RowColData rc, DofMethod &method, int lower_rank, int upper_rank, int verb=DEFAULT_VERBOSITY)=0
Make a loop over dofs.
virtual MoFEMErrorCode modify_problem_ref_level_add_bit(const std::string &name_problem, const BitRefLevel &bit)=0
add ref level to problem
DEPRECATED MoFEMErrorCode partition_problem(const std::string &name, int verb=-1)
partition problem dofs (collective)
DEPRECATED MoFEMErrorCode VecCreateGhost(const std::string &name, RowColData rc, Vec *V) const
create ghost vector for problem (collective)collective - need to be run on all processors in communic...
PetscErrorCode PetscOptionsGetInt(PetscOptions *, const char pre[], const char name[], PetscInt *ivalue, PetscBool *set)
#define MYPCOMM_INDEX
default communicator number PCOMM
Definition: definitions.h:292
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
DEPRECATED MoFEMErrorCode partition_ghost_dofs(const std::string &name, int verb=-1)
determine ghost nodes
virtual MoFEMErrorCode build_fields(int verb=DEFAULT_VERBOSITY)=0
MoFEMErrorCode generateReferenceElementMesh()
Generate reference mesh on single element.
virtual MoFEMErrorCode add_problem(const std::string &name, enum MoFEMTypes bh=MF_EXCL, int verb=DEFAULT_VERBOSITY)=0
Add problem.
virtual bool check_field(const std::string &name) const =0
check if field is in database
DeprecatedCoreInterface Interface
Definition: Interface.hpp:1943
Core (interface) class.
Definition: Core.hpp:77
continuous field
Definition: definitions.h:177
DEPRECATED MoFEMErrorCode MatCreateMPIAIJWithArrays(const std::string &name, Mat *Aij, int verb=DEFAULT_VERBOSITY)
DEPRECATED MoFEMErrorCode seed_ref_level_3D(const EntityHandle meshset, const BitRefLevel &bit, int verb=-1)
seed 2D entities in the meshset and their adjacencies (only TETs adjacencies) in a particular BitRefL...
DEPRECATED MoFEMErrorCode add_ents_to_field_by_TETs(const EntityHandle meshset, const std::string &name, int verb=-1)
set field entities from adjacencies of tetrahedronThe lower dimension entities are added depending on...
PetscErrorCode PetscOptionsGetBool(PetscOptions *, const char pre[], const char name[], PetscBool *bval, PetscBool *set)
virtual MoFEMErrorCode build_finite_elements(int verb=DEFAULT_VERBOSITY)=0
Build finite elementsBuild finite element data structures. Have to be run before problem and adjacenc...
DEPRECATED MoFEMErrorCode loop_finite_elements(const Problem *problem_ptr, const std::string &fe_name, FEMethod &method, int lower_rank, int upper_rank, MoFEMTypes bh, int verb=DEFAULT_VERBOSITY)