00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028 #ifndef _ABSTRACTNONLINEARASSEMBLER_HPP_
00029 #define _ABSTRACTNONLINEARASSEMBLER_HPP_
00030
00036 #include <vector>
00037
00038 #include "AbstractStaticAssembler.hpp"
00039 #include "BoundaryConditionsContainer.hpp"
00040 #include "AbstractNonlinearSolver.hpp"
00041 #include "AbstractNonlinearSolver.hpp"
00042 #include "SimplePetscNonlinearSolver.hpp"
00043 #include "RandomNumberGenerator.hpp"
00044 #include "PetscTools.hpp"
00045
00046
00047
00048
00049
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, class CONCRETE>
00062 PetscErrorCode AbstractNonlinearAssembler_AssembleResidual(SNES snes,
00063 Vec currentGuess,
00064 Vec residualVector,
00065 void *pContext);
00066
00067 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, class CONCRETE>
00068 PetscErrorCode AbstractNonlinearAssembler_AssembleJacobian(SNES snes,
00069 Vec currentGuess,
00070 Mat *pGlobalJacobian,
00071 Mat *pPreconditioner,
00072 MatStructure *pMatStructure,
00073 void *pContext);
00074
00075
00076
00080 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, class CONCRETE>
00081 class AbstractNonlinearAssembler : public AbstractStaticAssembler<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM,true,CONCRETE>
00082 {
00083 typedef AbstractStaticAssembler<ELEMENT_DIM,SPACE_DIM,PROBLEM_DIM,true,CONCRETE> BaseClassType;
00084 private:
00085 Vec mInitialGuess;
00086 protected:
00087 AbstractNonlinearSolver* mpSolver;
00088 bool mWeAllocatedSolverMemory;
00089
00090 bool mUseAnalyticalJacobian;
00091
00095 void ApplyDirichletConditions(Vec currentGuess, bool applyToMatrix)
00096 {
00097 Vec& residual = this->mpLinearSystem->rGetRhsVector();
00098 Mat& jacobian = this->mpLinearSystem->rGetLhsMatrix();
00099 assert((jacobian && applyToMatrix) || (!jacobian && !applyToMatrix));
00100
00101 if (residual)
00102 {
00103 this->mpBoundaryConditions->ApplyDirichletToNonlinearResidual(currentGuess, residual);
00104 }
00105 if (jacobian)
00106 {
00107 this->mpBoundaryConditions->ApplyDirichletToNonlinearJacobian(jacobian);
00108 }
00109 }
00110
00111 public :
00122 PetscErrorCode AssembleResidual(const Vec currentGuess, Vec residualVector)
00123 {
00124
00125
00126 this->PrepareForSolve();
00127 delete this->mpLinearSystem;
00128 this->mpLinearSystem = new LinearSystem(residualVector, NULL);
00129 this->AssembleSystem(true, false, currentGuess, 0.0);
00130 return 0;
00131 }
00132
00133
00146 PetscErrorCode AssembleJacobian(const Vec currentGuess, Mat *pGlobalJacobian)
00147 {
00148 if (mUseAnalyticalJacobian)
00149 {
00150
00151
00152 delete this->mpLinearSystem;
00153 this->mpLinearSystem = new LinearSystem(NULL, *pGlobalJacobian);
00154 this->AssembleSystem(false, true, currentGuess, 0.0);
00155 return 0;
00156 }
00157 else
00158 {
00159 return AssembleJacobianNumerically(currentGuess, pGlobalJacobian);
00160 }
00161 }
00162
00163 protected:
00173 PetscErrorCode AssembleJacobianNumerically(const Vec currentGuess, Mat *pJacobian)
00174 {
00175 unsigned num_nodes = PROBLEM_DIM*this->mpMesh->GetNumNodes();
00176
00177
00178 Vec residual;
00179 Vec perturbed_residual;
00180 Vec result;
00181
00182 VecCreate(PETSC_COMM_WORLD, &residual);
00183 VecCreate(PETSC_COMM_WORLD, &result);
00184 VecCreate(PETSC_COMM_WORLD, &perturbed_residual);
00185
00186 VecSetSizes(residual, PETSC_DECIDE,num_nodes);
00187 VecSetSizes(result, PETSC_DECIDE,num_nodes);
00188 VecSetSizes(perturbed_residual,PETSC_DECIDE,num_nodes);
00189
00190 VecSetFromOptions(residual);
00191 VecSetFromOptions(result);
00192 VecSetFromOptions(perturbed_residual);
00193
00194
00195 Vec current_guess_copy;
00196 PETSCEXCEPT( VecDuplicate(currentGuess, ¤t_guess_copy) );
00197 PETSCEXCEPT( VecCopy(currentGuess, current_guess_copy) );
00198
00199
00200 AssembleResidual(currentGuess, residual);
00201
00202
00203 double h = 0.00001;
00204 PetscScalar subtract = -1;
00205 PetscScalar one_over_h = 1.0/h;
00206
00207 PetscInt ilo, ihi;
00208 VecGetOwnershipRange(current_guess_copy, &ilo, &ihi);
00209 unsigned lo=ilo;
00210 unsigned hi=ihi;
00211
00212
00213 for (unsigned global_index_outer = 0; global_index_outer < num_nodes; global_index_outer++)
00214 {
00215
00216 if (lo<=global_index_outer && global_index_outer<hi)
00217 {
00218 PETSCEXCEPT( VecSetValue(current_guess_copy, global_index_outer,h, ADD_VALUES) );
00219 }
00220 AssembleResidual(current_guess_copy, perturbed_residual);
00221
00222
00223 #if (PETSC_VERSION_MINOR == 2) //Old API
00224 PETSCEXCEPT( VecWAXPY(&subtract, residual, perturbed_residual, result) );
00225 PETSCEXCEPT( VecScale(&one_over_h, result) );
00226 #else
00227 PETSCEXCEPT( VecWAXPY(result, subtract, residual, perturbed_residual) );
00228 PETSCEXCEPT( VecScale(result, one_over_h) );
00229 #endif
00230
00231 double *p_result;
00232 PETSCEXCEPT( VecGetArray(result, &p_result) );
00233 for (unsigned global_index=lo; global_index < hi; global_index++)
00234 {
00235 unsigned local_index = global_index - lo;
00236 PETSCEXCEPT( MatSetValue(*pJacobian, global_index, global_index_outer,
00237 p_result[local_index], INSERT_VALUES) );
00238 }
00239 PETSCEXCEPT( VecRestoreArray(result, &p_result) );
00240
00241 if (lo<=global_index_outer && global_index_outer<hi)
00242 {
00243 PETSCEXCEPT( VecSetValue(current_guess_copy, global_index_outer, -h, ADD_VALUES) );
00244 }
00245 }
00246
00247 VecDestroy(residual);
00248 VecDestroy(perturbed_residual);
00249 VecDestroy(result);
00250 VecDestroy(current_guess_copy);
00251
00252 MatAssemblyBegin(*pJacobian, MAT_FINAL_ASSEMBLY);
00253 MatAssemblyEnd(*pJacobian, MAT_FINAL_ASSEMBLY);
00254
00255 return 0;
00256 }
00257
00258 virtual void AssembleSystem(bool assembleVector, bool assembleMatrix,
00259 Vec currentGuess=NULL, double currentTime=0.0)
00260 {
00261
00262 assert( currentGuess );
00263 BaseClassType::AssembleSystem(assembleVector, assembleMatrix, currentGuess, currentTime);
00264 }
00265
00266 bool ProblemIsNonlinear()
00267 {
00268 return true;
00269 }
00270
00275 void InitialiseForSolve(Vec initialGuess)
00276 {
00277
00278 PetscInt size_of_init_guess;
00279 VecGetSize(initialGuess, &size_of_init_guess);
00280 PetscInt problem_size = PROBLEM_DIM * this->mpMesh->GetNumNodes();
00281 if (size_of_init_guess != problem_size)
00282 {
00283 std::stringstream error_message;
00284 error_message << "Size of initial guess vector, " << size_of_init_guess
00285 << ", does not match size of problem, " << problem_size;
00286 EXCEPTION(error_message.str());
00287 }
00288 }
00289
00298 Vec StaticSolve(Vec currentSolutionOrGuess=NULL,
00299 double currentTime=0.0,
00300 bool assembleMatrix=true)
00301 {
00302 assert(this->mpBoundaryConditions!=NULL);
00303 assert(currentSolutionOrGuess!=NULL);
00304 assert(assembleMatrix);
00305
00306
00307
00308 Vec answer = this->mpSolver->Solve(&AbstractNonlinearAssembler_AssembleResidual<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CONCRETE>,
00309 &AbstractNonlinearAssembler_AssembleJacobian<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CONCRETE>,
00310 currentSolutionOrGuess,
00311 this);
00312 return answer;
00313 }
00314
00315 public:
00319 AbstractNonlinearAssembler(unsigned numQuadPoints = 2) :
00320 BaseClassType(numQuadPoints)
00321 {
00322 mpSolver = new SimplePetscNonlinearSolver;
00323 mWeAllocatedSolverMemory = true;
00324
00325 mUseAnalyticalJacobian = false;
00326
00327 mInitialGuess = NULL;
00328 }
00329
00330 ~AbstractNonlinearAssembler()
00331 {
00332 if (mWeAllocatedSolverMemory)
00333 {
00334 delete mpSolver;
00335 }
00336 }
00337
00346 void SetUseAnalyticalJacobian(bool useAnalyticalJacobian)
00347 {
00348 mUseAnalyticalJacobian = useAnalyticalJacobian;
00349 }
00350
00359 virtual Vec Solve(Vec initialGuess, bool useAnalyticalJacobian=false)
00360 {
00361 assert(initialGuess!=NULL);
00362 SetUseAnalyticalJacobian(useAnalyticalJacobian);
00363
00364 this->PrepareForSolve();
00365 InitialiseForSolve(initialGuess);
00366
00367 return StaticSolve(initialGuess);
00368 }
00369
00375 void SetNonlinearSolver(AbstractNonlinearSolver* pNonlinearSolver)
00376 {
00377 if (mWeAllocatedSolverMemory)
00378 {
00379 delete mpSolver;
00380 }
00381 mpSolver = pNonlinearSolver;
00382 mWeAllocatedSolverMemory = false;
00383 }
00384
00388 Vec CreateConstantInitialGuess(double value)
00389 {
00390 assert(this->mpMesh!=NULL);
00391 unsigned size = PROBLEM_DIM * this->mpMesh->GetNumNodes();
00392 return PetscTools::CreateVec(size, value);
00393 }
00394
00411 bool VerifyJacobian(double tol=1e-4, bool print=false)
00412 {
00413 unsigned size = PROBLEM_DIM * this->mpMesh->GetNumNodes();
00414
00415 Vec initial_guess;
00416 VecCreate(PETSC_COMM_WORLD, &initial_guess);
00417 VecSetSizes(initial_guess, PETSC_DECIDE, size);
00418 VecSetFromOptions(initial_guess);
00419
00420
00421 for (unsigned i=0; i<size; i++)
00422 {
00423 VecSetValue(initial_guess, i, 0.0, INSERT_VALUES);
00424 }
00425
00426 VecAssemblyBegin(initial_guess);
00427 VecAssemblyEnd(initial_guess);
00428
00429
00430 Mat analytic_jacobian;
00431 Mat numerical_jacobian;
00432
00433 PetscTools::SetupMat(analytic_jacobian, size, size);
00434 PetscTools::SetupMat(numerical_jacobian, size, size);
00435
00436 mUseAnalyticalJacobian = true;
00437 AssembleJacobian(initial_guess, &analytic_jacobian);
00438
00439 mUseAnalyticalJacobian = false;
00440 AssembleJacobian(initial_guess, &numerical_jacobian);
00441
00442 bool all_less_than_tol = true;
00443
00444 if (print)
00445 {
00446 std::cout << "Difference between numerical and analyical Jacobians:\n\n";
00447 }
00448 for (unsigned i=0; i<size; i++)
00449 {
00450 for (unsigned j=0; j<size; j++)
00451 {
00452 double val_a[1];
00453 double val_n[1];
00454 PetscInt row[1];
00455 PetscInt col[1];
00456 row[0] = i;
00457 col[0] = j;
00458
00459 MatGetValues(numerical_jacobian,1,row,1,col,val_n);
00460 MatGetValues(analytic_jacobian,1,row,1,col,val_a);
00461
00462 if (print)
00463 {
00464 std::cout << val_n[0] - val_a[0]<< " ";
00465 }
00466
00467 if (fabs(val_n[0]-val_a[0]) > tol)
00468 {
00469 #define COVERAGE_IGNORE // would have to write a bad concrete assembler class just to cover this line
00470 all_less_than_tol = false;
00471 #undef COVERAGE_IGNORE
00472 }
00473 }
00474 if (print)
00475 {
00476 std::cout << "\n" <<std::flush;
00477 }
00478 }
00479 std::cout << std::flush;
00480
00481 MatDestroy(numerical_jacobian);
00482 MatDestroy(analytic_jacobian);
00483 VecDestroy(initial_guess);
00484
00485 return all_less_than_tol;
00486 }
00487
00488 };
00489
00490
00491
00492
00493
00494
00495
00496
00497
00498
00499
00500
00501
00502
00503
00504
00518 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, class CONCRETE>
00519 PetscErrorCode AbstractNonlinearAssembler_AssembleResidual(SNES snes, Vec currentGuess,
00520 Vec residualVector, void *pContext)
00521 {
00522
00523 AbstractNonlinearAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CONCRETE> *pAssembler =
00524 (AbstractNonlinearAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CONCRETE>*) pContext;
00525
00526 PetscErrorCode ierr = pAssembler->AssembleResidual(currentGuess, residualVector);
00527
00528
00529
00530
00531
00532 return ierr;
00533 }
00534
00535
00536
00552 template<unsigned ELEMENT_DIM, unsigned SPACE_DIM, unsigned PROBLEM_DIM, class CONCRETE>
00553 PetscErrorCode AbstractNonlinearAssembler_AssembleJacobian(SNES snes, Vec currentGuess,
00554 Mat *pGlobalJacobian, Mat *pPreconditioner,
00555 MatStructure *pMatStructure, void *pContext)
00556 {
00557
00558
00559
00560 AbstractNonlinearAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CONCRETE> *pAssembler =
00561 (AbstractNonlinearAssembler<ELEMENT_DIM, SPACE_DIM, PROBLEM_DIM, CONCRETE>*) pContext;
00562
00563 PetscErrorCode ierr = pAssembler->AssembleJacobian(currentGuess, pGlobalJacobian);
00564
00565
00566 return ierr;
00567 }
00568
00569
00570
00571 #endif //_ABSTRACTNONLINEARASSEMBLER_HPP_