v0.16.3
Loading...
Searching...
No Matches
Public Member Functions | Private Attributes | Friends | List of all members
EshelbianPlasticity::TopologicalTAOCtxImpl Struct Reference
Inheritance diagram for EshelbianPlasticity::TopologicalTAOCtxImpl:
[legend]
Collaboration diagram for EshelbianPlasticity::TopologicalTAOCtxImpl:
[legend]

Public Member Functions

 TopologicalTAOCtxImpl (EshelbianCore *ep, ForcesAndSourcesCore::GaussHookFun set_integration_at_interior, ForcesAndSourcesCore::GaussHookFun set_integration_at_face, SmartPetscObj< TS > time_solver)
 
- Public Member Functions inherited from EshelbianPlasticity::TopologicalTAOCtx
 TopologicalTAOCtx ()=default
 
virtual ~TopologicalTAOCtx ()=default
 

Private Attributes

EshelbianCore * ep_ptr
 
ForcesAndSourcesCore::GaussHookFun integrationAtInterior
 
ForcesAndSourcesCore::GaussHookFun integrationAtFace
 
SmartPetscObj< TS > timeSolver
 
SmartPetscObj< Vec > primalProblemVec
 
SmartPetscObj< Vec > adjointProblemVec
 

Friends

MoFEMErrorCode solveEquilibriumStateTS (TopologicalTAOCtxImpl *ctx_impl_ptr)
 
MoFEMErrorCode solveAdjointSensitivityKSP (TopologicalTAOCtxImpl *ctx_impl_ptr, ObjectiveModelType eval_energy_model)
 
MoFEMErrorCode evaluateGradientImpl (TopologicalTAOCtxImpl *ctx_impl_ptr, double *f, Vec g, ObjectiveModelType eval_energy_model)
 
MoFEMErrorCode evaluateObjectiveImpl (TopologicalTAOCtxImpl *ctx_impl_ptr, double *f, ObjectiveModelType eval_energy_model)
 
MoFEMErrorCode finiteDifferenceGradientTest (TopologicalTAOCtxImpl *ctx_impl_ptr, Vec sol, double *f, Vec g, double epsilon, ObjectiveModelType eval_energy_model)
 
MoFEMErrorCode finiteDifference_dJdX_Test (TopologicalTAOCtxImpl *ctx_impl_ptr, double epsilon, ObjectiveModelType eval_energy_model)
 
MoFEMErrorCode finiteDifference_dJdx_Test (TopologicalTAOCtxImpl *ctx_impl_ptr, double epsilon, ObjectiveModelType eval_energy_model)
 
MoFEMErrorCode finiteDifference_dJd_adjoint_Test (TopologicalTAOCtxImpl *ctx_impl_ptr, double epsilon)
 
MoFEMErrorCode topologicalEvaluateObjectiveAndGradient (Tao tao, Vec sol, PetscReal *f, Vec g, void *ctx)
 
MoFEMErrorCode topologicalEvaluateObjectiveAndGradient (Tao tao, Vec sol, PetscReal *f, Vec g, void *ctx)
 
MoFEMErrorCode testTopologicalDerivative (TopologicalTAOCtx *ctx_ptr, Vec sol, PetscReal *f, Vec g, ObjectiveModelType eval_energy_model)
 

Detailed Description

Definition at line 32 of file EshelbianTopologicalDerivative.cpp.

Constructor & Destructor Documentation

◆ TopologicalTAOCtxImpl()

EshelbianPlasticity::TopologicalTAOCtxImpl::TopologicalTAOCtxImpl ( EshelbianCore *  ep,
ForcesAndSourcesCore::GaussHookFun  set_integration_at_interior,
ForcesAndSourcesCore::GaussHookFun  set_integration_at_face,
SmartPetscObj< TS >  time_solver 
)
inline

Definition at line 33 of file EshelbianTopologicalDerivative.cpp.

38 : ep_ptr(ep), integrationAtInterior(set_integration_at_interior),
39 integrationAtFace(set_integration_at_face), timeSolver(time_solver) {
40
41 primalProblemVec = createDMVector(ep->dmElastic, RowColData::COL);
42 adjointProblemVec = createDMVector(ep->dmElastic, RowColData::COL);
43 }
COL
@ COL
Definition definitions.h:136
MoFEM::createDMVector
auto createDMVector(DM dm, RowColData rc=RowColData::COL)
Get smart vector from DM.
Definition DMMoFEM.hpp:1237
EshelbianCore::dmElastic
SmartPetscObj< DM > dmElastic
Elastic problem.
Definition EshelbianCore.hpp:232

Friends And Related Symbol Documentation

◆ evaluateGradientImpl

MoFEMErrorCode evaluateGradientImpl ( TopologicalTAOCtxImpl *  ctx_impl_ptr,
double *  f,
Vec  g,
ObjectiveModelType  eval_energy_model 
)
friend

Definition at line 201 of file EshelbianTopologicalDerivative.cpp.

203 {
204 MoFEMFunctionBegin;
205 auto &ep = *ctx_impl_ptr->ep_ptr;
206 auto &m_field = ep.mField;
207 auto interior_integration_hook = ctx_impl_ptr->integrationAtInterior;
208 auto boundary_integration_hook = ctx_impl_ptr->integrationAtFace;
209
210 auto fe_material =
211 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
212 auto fe_interior_adjoint =
213 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
214 auto fe_natural_adjoint =
215 boost::make_shared<FaceElementForcesAndSourcesCore>(m_field);
216
217 boost::shared_ptr<double> J_ptr(
218 f, [](double *) {}); // Custom deleter to avoid deallocation
219 auto dJ_dX_vec = SmartPetscObj<Vec>(g, true);
220 CHKERR VecZeroEntries(dJ_dX_vec);
221 CHKERR VecGhostUpdateBegin(dJ_dX_vec, INSERT_VALUES, SCATTER_FORWARD);
222 CHKERR VecGhostUpdateEnd(dJ_dX_vec, INSERT_VALUES, SCATTER_FORWARD);
223
224 CHKERR pushTopologicalMaterialOps(ep, fe_material, interior_integration_hook,
225 boundary_integration_hook, J_ptr,
226 dJ_dX_vec, eval_energy_model);
227
228 CHKERR VecScale(ctx_impl_ptr->adjointProblemVec, -1.0);
229 auto alpha = ep.alphaW;
230 auto rho = ep.alphaRho;
231 auto alpha_viscous_omega = ep.alphaViscousOmega;
232 CHKERR pushTopologicalInteriorOps_dJ_adjoint_gradient(
233 ep, fe_interior_adjoint, interior_integration_hook,
234 boundary_integration_hook, ctx_impl_ptr->adjointProblemVec, dJ_dX_vec,
235 alpha, rho, alpha_viscous_omega, nullptr);
236 CHKERR pushTopologicalBoundaryOps_dJ_adjoint_gradient(
237 ep, fe_natural_adjoint, interior_integration_hook,
238 boundary_integration_hook, ctx_impl_ptr->adjointProblemVec, dJ_dX_vec,
239 nullptr);
240
241 CHKERR DMoFEMLoopFiniteElements(ep.dmMaterial, ep.elementVolumeName,
242 fe_material);
243 CHKERR DMoFEMLoopFiniteElements(ep.dmMaterial, ep.elementVolumeName,
244 fe_interior_adjoint);
245 CHKERR DMoFEMLoopFiniteElements(ep.dmMaterial, ep.naturalBcElement,
246 fe_natural_adjoint);
247
248 CHKERR VecAssemblyBegin(dJ_dX_vec);
249 CHKERR VecAssemblyEnd(dJ_dX_vec);
250 CHKERR VecGhostUpdateBegin(dJ_dX_vec, ADD_VALUES, SCATTER_REVERSE);
251 CHKERR VecGhostUpdateEnd(dJ_dX_vec, ADD_VALUES, SCATTER_REVERSE);
252 CHKERR VecGhostUpdateBegin(dJ_dX_vec, INSERT_VALUES, SCATTER_FORWARD);
253 CHKERR VecGhostUpdateEnd(dJ_dX_vec, INSERT_VALUES, SCATTER_FORWARD);
254
255 MoFEMFunctionReturn(0);
256}
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
MoFEM::DMoFEMLoopFiniteElements
PetscErrorCode DMoFEMLoopFiniteElements(DM dm, const char fe_name[], MoFEM::FEMethod *method, CacheTupleWeakPtr cache_ptr=CacheTupleSharedPtr())
Executes FEMethod for finite elements in DM.
Definition DMMoFEM.cpp:576
EshelbianPlasticity::pushTopologicalBoundaryOps_dJ_adjoint_gradient
MoFEMErrorCode pushTopologicalBoundaryOps_dJ_adjoint_gradient(EshelbianCore &ep, boost::shared_ptr< FaceElementForcesAndSourcesCore > fe, ForcesAndSourcesCore::GaussHookFun interior_integration_hook, ForcesAndSourcesCore::GaussHookFun boundary_integration_hook, SmartPetscObj< Vec > lambda_vec, SmartPetscObj< Vec > dJ_dX_vec, boost::shared_ptr< double > J_ptr, SmartPetscObj< Vec > topo_vec=SmartPetscObj< Vec >())
Definition EshelbianTopologicalDerivativePush.cpp:464
EshelbianPlasticity::pushTopologicalInteriorOps_dJ_adjoint_gradient
MoFEMErrorCode pushTopologicalInteriorOps_dJ_adjoint_gradient(EshelbianCore &ep, boost::shared_ptr< VolumeElementForcesAndSourcesCore > fe, ForcesAndSourcesCore::GaussHookFun interior_integration_hook, ForcesAndSourcesCore::GaussHookFun boundary_integration_hook, SmartPetscObj< Vec > lambda_vec, SmartPetscObj< Vec > dJ_dX_vec, const double alpha, const double rho, const double alpha_viscous_omega, boost::shared_ptr< double > J_ptr, SmartPetscObj< Vec > topo_vec=SmartPetscObj< Vec >())
Definition EshelbianTopologicalDerivativePush.cpp:306
EshelbianPlasticity::pushTopologicalMaterialOps
MoFEMErrorCode pushTopologicalMaterialOps(EshelbianCore &ep, boost::shared_ptr< VolumeElementForcesAndSourcesCore > fe, ForcesAndSourcesCore::GaussHookFun interior_integration_hook, ForcesAndSourcesCore::GaussHookFun boundary_integration_hook, boost::shared_ptr< double > J_ptr, SmartPetscObj< Vec > dJ_dX_vec, ObjectiveModelType eval_energy_model)
Definition EshelbianTopologicalDerivativePush.cpp:251
g
constexpr double g
Definition shallow_wave.cpp:63
EshelbianCore::mField
MoFEM::Interface & mField
Definition EshelbianCore.hpp:218
MoFEM::SmartPetscObj
intrusive_ptr for managing petsc objects
Definition PetscSmartObj.hpp:85
rho
double rho
Definition plastic.cpp:144

◆ evaluateObjectiveImpl

MoFEMErrorCode evaluateObjectiveImpl ( TopologicalTAOCtxImpl *  ctx_impl_ptr,
double *  f,
ObjectiveModelType  eval_energy_model 
)
friend

Definition at line 258 of file EshelbianTopologicalDerivative.cpp.

260 {
261 MoFEMFunctionBegin;
262 auto &ep = *ctx_impl_ptr->ep_ptr;
263 auto &m_field = ep.mField;
264 auto interior_integration_hook = ctx_impl_ptr->integrationAtInterior;
265 auto boundary_integration_hook = ctx_impl_ptr->integrationAtFace;
266
267 auto fe_material =
268 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
269 *f = 0;
270 boost::shared_ptr<double> J_ptr(
271 f, [](double *) {}); // Custom deleter to avoid deallocation
272 CHKERR pushTopologicalMaterialOps(ep, fe_material, interior_integration_hook,
273 boundary_integration_hook, J_ptr,
274 SmartPetscObj<Vec>(), eval_energy_model);
275 CHKERR DMoFEMLoopFiniteElements(ep.dmMaterial, ep.elementVolumeName,
276 fe_material);
277
278 MoFEMFunctionReturn(0);
279}

◆ finiteDifference_dJd_adjoint_Test

MoFEMErrorCode finiteDifference_dJd_adjoint_Test ( TopologicalTAOCtxImpl *  ctx_impl_ptr,
double  epsilon 
)
friend

Definition at line 768 of file EshelbianTopologicalDerivative.cpp.

769 {
770 MoFEMFunctionBegin;
771 auto &ep = *ctx_impl_ptr->ep_ptr;
772
773 auto row_sol =
774 createDMVector(ctx_impl_ptr->ep_ptr->dmMaterial, RowColData::ROW);
775 CHKERR DMoFEMMeshToLocalVector(ep.dmMaterial, row_sol, INSERT_VALUES,
776 SCATTER_FORWARD, RowColData::ROW);
777 CHKERR VecGhostUpdateBegin(row_sol, INSERT_VALUES, SCATTER_FORWARD);
778 CHKERR VecGhostUpdateEnd(row_sol, INSERT_VALUES, SCATTER_FORWARD);
779
780 auto col_sol =
781 createDMVector(ctx_impl_ptr->ep_ptr->dmElastic, RowColData::COL);
782 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, col_sol, INSERT_VALUES,
783 SCATTER_FORWARD, RowColData::COL);
784 CHKERR VecGhostUpdateBegin(col_sol, INSERT_VALUES, SCATTER_FORWARD);
785 CHKERR VecGhostUpdateEnd(col_sol, INSERT_VALUES, SCATTER_FORWARD);
786
787 MOFEM_LOG("EP", Sev::inform) << "Starting finite difference dJ_adjoint_dX "
788 "gradient test with epsilon = "
789 << epsilon;
790
791 auto opt = ep.mField.getInterface<OperatorsTester>();
792 auto direction_vec =
793 opt->setRandomFields(ep.dmMaterial, {{ep.materialH1Positions, {-1., 1.}}},
794 nullptr, RowColData::ROW);
795 CHKERR VecGhostUpdateBegin(direction_vec, INSERT_VALUES, SCATTER_FORWARD);
796 CHKERR VecGhostUpdateEnd(direction_vec, INSERT_VALUES, SCATTER_FORWARD);
797
798 const std::array<double, 2> piola_range{{-1, 1}};
799 const std::array<double, 2> bubble_range{{-1, 1}};
800 const std::array<double, 2> spatial_l2_disp_range{{-1, 1}};
801 const std::array<double, 2> rot_axis_range{{-1, 1}};
802 const std::array<double, 2> stretch_tensor_range{{-1, 1}};
803 const std::array<double, 2> hybrid_spatial_disp_range{{-1, 1}};
804
805 std::vector<OperatorsTester::RandomFieldData> adjoint_random_fields{
806 {ep.piolaStress, piola_range},
807 {ep.bubbleField, bubble_range},
808 {ep.spatialL2Disp, spatial_l2_disp_range},
809 {ep.rotAxis, rot_axis_range},
810 {ep.stretchTensor, stretch_tensor_range},
811 {ep.hybridSpatialDisp, hybrid_spatial_disp_range}};
812
813 auto adjoint_vec = opt->setRandomFields(ep.dmMaterial, adjoint_random_fields,
814 nullptr, RowColData::COL);
815 CHKERR VecGhostUpdateBegin(adjoint_vec, INSERT_VALUES, SCATTER_FORWARD);
816 CHKERR VecGhostUpdateEnd(adjoint_vec, INSERT_VALUES, SCATTER_FORWARD);
817
818 const std::array<double, 2> state_piola_range{{-1, 1}};
819 const std::array<double, 2> state_bubble_range{{-1, 1}};
820 const std::array<double, 2> state_spatial_l2_disp_range{{-1, 1}};
821 const std::array<double, 2> state_rot_axis_range{{-1, 1}};
822 const std::array<double, 2> state_stretch_tensor_range{{-1, 1}};
823 const std::array<double, 2> state_hybrid_spatial_disp_range{{-1, 1}};
824
825 std::vector<OperatorsTester::RandomFieldData> state_random_fields{
826 {ep.piolaStress, state_piola_range},
827 {ep.bubbleField, state_bubble_range},
828 {ep.spatialL2Disp, state_spatial_l2_disp_range},
829 {ep.rotAxis, state_rot_axis_range},
830 {ep.stretchTensor, state_stretch_tensor_range},
831 {ep.hybridSpatialDisp, state_hybrid_spatial_disp_range}};
832
833 auto state_vec = opt->setRandomFields(ep.dmElastic, state_random_fields,
834 nullptr, RowColData::COL);
835 CHKERR VecGhostUpdateBegin(state_vec, INSERT_VALUES, SCATTER_FORWARD);
836 CHKERR VecGhostUpdateEnd(state_vec, INSERT_VALUES, SCATTER_FORWARD);
837
838 CHKERR setSpatialConfigurationImpl(ep, state_vec);
839
840 boost::shared_ptr<double> J_ptr;
841 auto dJ_dX_vec = vectorDuplicate(direction_vec);
842 CHKERR VecZeroEntries(dJ_dX_vec);
843 auto x_t_vec = vectorDuplicate(adjoint_vec);
844 CHKERR VecZeroEntries(x_t_vec);
845
846 auto get_f_rhs_vec = [&]() {
847 auto f_rhs_vec = vectorDuplicate(adjoint_vec);
848 auto fe_interior_rhs = ep.elasticFeRhs;
849 auto fe_boundary_rhs = ep.elasticBcRhs;
850 fe_interior_rhs->f = f_rhs_vec;
851 fe_interior_rhs->x_t = x_t_vec;
852 fe_interior_rhs->data_ctx |=
853 PetscData::CTX_SET_X | PetscData::CTX_SET_X_T | PetscData::CTX_SET_TIME;
854 CHKERR TSGetTime(ctx_impl_ptr->timeSolver, &fe_interior_rhs->ts_t);
855 CHKERR TSGetTimeStep(ctx_impl_ptr->timeSolver, &fe_interior_rhs->ts_dt);
856 CHKERR TSGetStepNumber(ctx_impl_ptr->timeSolver, &fe_interior_rhs->ts_step);
857 fe_boundary_rhs->f = f_rhs_vec;
858 fe_boundary_rhs->x_t = x_t_vec;
859 fe_boundary_rhs->data_ctx |=
860 PetscData::CTX_SET_X | PetscData::CTX_SET_X_T | PetscData::CTX_SET_TIME;
861 CHKERR TSGetTime(ctx_impl_ptr->timeSolver, &fe_boundary_rhs->ts_t);
862 CHKERR TSGetTimeStep(ctx_impl_ptr->timeSolver, &fe_boundary_rhs->ts_dt);
863 CHKERR TSGetStepNumber(ctx_impl_ptr->timeSolver, &fe_boundary_rhs->ts_step);
864 CHKERR VecZeroEntries(f_rhs_vec);
865 CHKERR DMoFEMLoopFiniteElements(ep.dmElastic, ep.elementVolumeName,
866 fe_interior_rhs);
867 CHKERR DMoFEMLoopFiniteElements(ep.dmElastic, ep.naturalBcElement,
868 fe_boundary_rhs);
869 CHKERR VecAssemblyBegin(f_rhs_vec);
870 CHKERR VecAssemblyEnd(f_rhs_vec);
871 CHKERR VecGhostUpdateBegin(f_rhs_vec, ADD_VALUES, SCATTER_REVERSE);
872 CHKERR VecGhostUpdateEnd(f_rhs_vec, ADD_VALUES, SCATTER_REVERSE);
873 return f_rhs_vec;
874 };
875
876 auto alpha = ep.alphaW;
877 auto rho = ep.alphaRho;
878 auto alpha_viscous_omega = ep.alphaViscousOmega;
879
880 auto interior_integration_hook = ctx_impl_ptr->integrationAtInterior;
881 auto boundary_integration_hook = ctx_impl_ptr->integrationAtFace;
882 auto fe_adjoint =
883 boost::make_shared<VolumeElementForcesAndSourcesCore>(ep.mField);
884 auto fe_adjoint_boundary =
885 boost::make_shared<FaceElementForcesAndSourcesCore>(ep.mField);
886 CHKERR pushTopologicalInteriorOps_dJ_adjoint_gradient(
887 ep, fe_adjoint, interior_integration_hook, boundary_integration_hook,
888 adjoint_vec, dJ_dX_vec, alpha, rho, alpha_viscous_omega,
889 J_ptr);
890 CHKERR pushTopologicalBoundaryOps_dJ_adjoint_gradient(
891 ep, fe_adjoint_boundary, interior_integration_hook,
892 boundary_integration_hook, adjoint_vec, dJ_dX_vec, J_ptr);
893
894 CHKERR DMoFEMLoopFiniteElements(ep.dmMaterial, ep.elementVolumeName,
895 fe_adjoint);
896 CHKERR DMoFEMLoopFiniteElements(ep.dmMaterial, ep.naturalBcElement,
897 fe_adjoint_boundary);
898 CHKERR VecAssemblyBegin(dJ_dX_vec);
899 CHKERR VecAssemblyEnd(dJ_dX_vec);
900 CHKERR VecGhostUpdateBegin(dJ_dX_vec, ADD_VALUES, SCATTER_REVERSE);
901 CHKERR VecGhostUpdateEnd(dJ_dX_vec, ADD_VALUES, SCATTER_REVERSE);
902 CHKERR VecGhostUpdateBegin(dJ_dX_vec, INSERT_VALUES, SCATTER_FORWARD);
903 CHKERR VecGhostUpdateEnd(dJ_dX_vec, INSERT_VALUES, SCATTER_FORWARD);
904
905 auto delta_vec = vectorDuplicate(direction_vec);
906 CHKERR VecCopy(direction_vec, delta_vec);
907 CHKERR VecScale(delta_vec, epsilon);
908
909 auto a_vec = vectorDuplicate(delta_vec);
910 CHKERR VecCopy(row_sol, a_vec);
911 CHKERR VecAXPY(a_vec, 1, delta_vec);
912 auto b_vec = vectorDuplicate(delta_vec);
913 CHKERR VecCopy(row_sol, b_vec);
914 CHKERR VecAXPY(b_vec, -1, delta_vec);
915
916 auto fe_adjoint_fd_plus =
917 boost::make_shared<VolumeElementForcesAndSourcesCore>(ep.mField);
918 auto fe_adjoint_fd_minus =
919 boost::make_shared<VolumeElementForcesAndSourcesCore>(ep.mField);
920 auto fe_adjoint_fd_plus_boundary =
921 boost::make_shared<FaceElementForcesAndSourcesCore>(ep.mField);
922 auto fe_adjoint_fd_minus_boundary =
923 boost::make_shared<FaceElementForcesAndSourcesCore>(ep.mField);
924
925 auto J_plus = boost::make_shared<double>(0);
926 CHKERR pushTopologicalInteriorOps_dJ_adjoint_gradient(
927 ep, fe_adjoint_fd_plus, interior_integration_hook,
928 boundary_integration_hook, adjoint_vec,
929 SmartPetscObj<Vec>(), alpha, rho, alpha_viscous_omega, J_plus, a_vec);
930 CHKERR pushTopologicalBoundaryOps_dJ_adjoint_gradient(
931 ep, fe_adjoint_fd_plus_boundary, interior_integration_hook,
932 boundary_integration_hook, adjoint_vec, SmartPetscObj<Vec>(), J_plus,
933 a_vec);
934 auto J_minus = boost::make_shared<double>(0);
935 CHKERR pushTopologicalInteriorOps_dJ_adjoint_gradient(
936 ep, fe_adjoint_fd_minus, interior_integration_hook,
937 boundary_integration_hook, adjoint_vec,
938 SmartPetscObj<Vec>(), alpha, rho, alpha_viscous_omega, J_minus, b_vec);
939 CHKERR pushTopologicalBoundaryOps_dJ_adjoint_gradient(
940 ep, fe_adjoint_fd_minus_boundary, interior_integration_hook,
941 boundary_integration_hook, adjoint_vec, SmartPetscObj<Vec>(), J_minus,
942 b_vec);
943
944 CHKERR setMaterialConfigurationImpl(ep, a_vec);
945 CHKERR DMoFEMLoopFiniteElements(ep.dmMaterial, ep.elementVolumeName,
946 fe_adjoint_fd_plus);
947 CHKERR DMoFEMLoopFiniteElements(ep.dmMaterial, ep.naturalBcElement,
948 fe_adjoint_fd_plus_boundary);
949 auto f_rhs_plus_vec = get_f_rhs_vec();
950 double nrm_f_rhs_plus;
951 CHKERR VecNorm(f_rhs_plus_vec, NORM_2, &nrm_f_rhs_plus);
952
953 MOFEM_LOG("EP", Sev::inform)
954 << "Objective at a_vec: J_plus = " << std::setprecision(12) << *J_plus
955 << ", Norm of f_rhs_plus_vec = " << nrm_f_rhs_plus;
956
957 CHKERR setMaterialConfigurationImpl(ep, b_vec);
958 CHKERR DMoFEMLoopFiniteElements(ep.dmMaterial, ep.elementVolumeName,
959 fe_adjoint_fd_minus);
960 CHKERR DMoFEMLoopFiniteElements(ep.dmMaterial, ep.naturalBcElement,
961 fe_adjoint_fd_minus_boundary);
962 auto f_rhs_minus_vec = get_f_rhs_vec();
963 double nrm_f_rhs_minus;
964 CHKERR VecNorm(f_rhs_minus_vec, NORM_2, &nrm_f_rhs_minus);
965
966 MOFEM_LOG("EP", Sev::inform)
967 << "Objective at b_vec: J_minus = " << std::setprecision(12) << *J_minus
968 << ", Norm of f_rhs_minus_vec = " << nrm_f_rhs_minus;
969
970 double dJ_da = (*J_plus - *J_minus) / (2 * epsilon);
971 CHKERR VecAXPY(f_rhs_plus_vec, -1, f_rhs_minus_vec);
972 CHKERR VecScale(f_rhs_plus_vec, 1. / (2 * epsilon));
973 double lambda_dJ_dX;
974 CHKERR VecDot(f_rhs_plus_vec, adjoint_vec, &lambda_dJ_dX);
975
976 double exact_dJ = 0;
977 CHKERR VecDot(dJ_dX_vec, direction_vec, &exact_dJ);
978 double fd_error = dJ_da - exact_dJ;
979 double fd_dJ_dX_error = lambda_dJ_dX - exact_dJ;
980
981 MOFEM_LOG("EP", Sev::inform)
982 << "Fd dJ_adjoint/dX = " << dJ_da << ", exact dJ/dX = " << exact_dJ
983 << ", error = " << fd_error;
984 MOFEM_LOG("EP", Sev::inform)
985 << "Fd lambda * dJ/dX = " << lambda_dJ_dX
986 << ", exact dJ/dX = " << exact_dJ << ", error = " << fd_dJ_dX_error;
987
988 CHKERR setMaterialConfigurationImpl(ep, row_sol);
989 CHKERR setSpatialConfigurationImpl(ep, col_sol);
990
991 MoFEMFunctionReturn(0);
992}
ROW
@ ROW
Definition definitions.h:136
MoFEM::DMoFEMMeshToLocalVector
PetscErrorCode DMoFEMMeshToLocalVector(DM dm, Vec l, InsertMode mode, ScatterMode scatter_mode, RowColData rc=RowColData::COL)
set local (or ghosted) vector values on mesh for partition only
Definition DMMoFEM.cpp:514
MOFEM_LOG
#define MOFEM_LOG(channel, severity)
Log.
Definition LogManager.hpp:316
EshelbianPlasticity::setSpatialConfigurationImpl
static MoFEMErrorCode setSpatialConfigurationImpl(EshelbianCore &ep, Vec sol)
Definition EshelbianTopologicalDerivative.cpp:290
EshelbianPlasticity::setMaterialConfigurationImpl
static MoFEMErrorCode setMaterialConfigurationImpl(EshelbianCore &ep, Vec sol)
Definition EshelbianTopologicalDerivative.cpp:281
MoFEM::vectorDuplicate
SmartPetscObj< Vec > vectorDuplicate(Vec vec)
Create duplicate vector of smart vector.
Definition PetscSmartObj.hpp:223
EshelbianCore::dmMaterial
SmartPetscObj< DM > dmMaterial
Material problem.
Definition EshelbianCore.hpp:233
MoFEM::OperatorsTester
Calculate directional derivative of the right hand side and compare it with tangent matrix derivative...
Definition OperatorsTester.hpp:21
MoFEM::OperatorsTester::setRandomFields
SmartPetscObj< Vec > setRandomFields(SmartPetscObj< DM > dm, std::vector< RandomFieldData > random_fields, boost::shared_ptr< Range > ents=nullptr, RowColData r=RowColData::COL)
Generate random fields.
Definition OperatorsTester.cpp:23
MoFEM::PetscData::CTX_SET_X_T
@ CTX_SET_X_T
Time derivative X_t is set.
Definition LoopMethods.hpp:94
MoFEM::PetscData::CTX_SET_TIME
@ CTX_SET_TIME
Time value is set.
Definition LoopMethods.hpp:96
MoFEM::PetscData::CTX_SET_X
@ CTX_SET_X
Solution vector X is set.
Definition LoopMethods.hpp:92

◆ finiteDifference_dJdX_Test

MoFEMErrorCode finiteDifference_dJdX_Test ( TopologicalTAOCtxImpl *  ctx_impl_ptr,
double  epsilon,
ObjectiveModelType  eval_energy_model 
)
friend

Definition at line 616 of file EshelbianTopologicalDerivative.cpp.

618 {
619 MoFEMFunctionBegin;
620 auto &ep = *ctx_impl_ptr->ep_ptr;
621 auto interior_integration_hook = ctx_impl_ptr->integrationAtInterior;
622 auto boundary_integration_hook = ctx_impl_ptr->integrationAtFace;
623
624 auto row_sol =
625 createDMVector(ctx_impl_ptr->ep_ptr->dmMaterial, RowColData::ROW);
626 CHKERR DMoFEMMeshToLocalVector(ep.dmMaterial, row_sol, INSERT_VALUES,
627 SCATTER_FORWARD, RowColData::ROW);
628 CHKERR VecGhostUpdateBegin(row_sol, INSERT_VALUES, SCATTER_FORWARD);
629 CHKERR VecGhostUpdateEnd(row_sol, INSERT_VALUES, SCATTER_FORWARD);
630
631 auto col_sol =
632 createDMVector(ctx_impl_ptr->ep_ptr->dmElastic, RowColData::COL);
633 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, col_sol, INSERT_VALUES,
634 SCATTER_FORWARD, RowColData::COL);
635 CHKERR VecGhostUpdateBegin(col_sol, INSERT_VALUES, SCATTER_FORWARD);
636 CHKERR VecGhostUpdateEnd(col_sol, INSERT_VALUES, SCATTER_FORWARD);
637
638 MOFEM_LOG("EP", Sev::inform) << "Starting finite difference dJ_dX "
639 "gradient test with epsilon = "
640 << epsilon;
641
642 auto opt = ep.mField.getInterface<OperatorsTester>();
643 auto direction_vec =
644 opt->setRandomFields(ep.dmMaterial, {{ep.materialH1Positions, {-1., 1.}}},
645 nullptr, RowColData::ROW);
646 CHKERR VecGhostUpdateBegin(direction_vec, INSERT_VALUES, SCATTER_FORWARD);
647 CHKERR VecGhostUpdateEnd(direction_vec, INSERT_VALUES, SCATTER_FORWARD);
648
649 const std::array<double, 2> piola_range{{-1, 1}};
650 const std::array<double, 2> bubble_range{{-1, 1}};
651 const std::array<double, 2> spatial_l2_disp_range{{-1, 1}};
652 const std::array<double, 2> rot_axis_range{{-1, 1}};
653 const std::array<double, 2> stretch_tensor_range{{-1, 1}};
654 const std::array<double, 2> hybrid_spatial_disp_range{{-1, 1}};
655
656 std::vector<OperatorsTester::RandomFieldData> adjoint_random_fields{
657 {ep.piolaStress, piola_range},
658 {ep.bubbleField, bubble_range},
659 {ep.spatialL2Disp, spatial_l2_disp_range},
660 {ep.rotAxis, rot_axis_range},
661 {ep.stretchTensor, stretch_tensor_range},
662 {ep.hybridSpatialDisp, hybrid_spatial_disp_range}};
663
664 auto adjoint_vec = opt->setRandomFields(ep.dmMaterial, adjoint_random_fields,
665 nullptr, RowColData::COL);
666 CHKERR VecGhostUpdateBegin(adjoint_vec, INSERT_VALUES, SCATTER_FORWARD);
667 CHKERR VecGhostUpdateEnd(adjoint_vec, INSERT_VALUES, SCATTER_FORWARD);
668
669 const std::array<double, 2> state_piola_range{{-1, 1}};
670 const std::array<double, 2> state_bubble_range{{-1, 1}};
671 const std::array<double, 2> state_spatial_l2_disp_range{{-1, 1}};
672 const std::array<double, 2> state_rot_axis_range{{-1, 1}};
673 const std::array<double, 2> state_stretch_tensor_range{{-1, 1}};
674 const std::array<double, 2> state_hybrid_spatial_disp_range{{-1, 1}};
675
676 std::vector<OperatorsTester::RandomFieldData> state_random_fields{
677 {ep.piolaStress, state_piola_range},
678 {ep.bubbleField, state_bubble_range},
679 {ep.spatialL2Disp, state_spatial_l2_disp_range},
680 {ep.rotAxis, state_rot_axis_range},
681 {ep.stretchTensor, state_stretch_tensor_range},
682 {ep.hybridSpatialDisp, state_hybrid_spatial_disp_range}};
683
684 auto state_vec = opt->setRandomFields(ep.dmElastic, state_random_fields,
685 nullptr, RowColData::COL);
686 CHKERR VecGhostUpdateBegin(state_vec, INSERT_VALUES, SCATTER_FORWARD);
687 CHKERR VecGhostUpdateEnd(state_vec, INSERT_VALUES, SCATTER_FORWARD);
688
689 CHKERR setSpatialConfigurationImpl(ep, state_vec);
690
691 boost::shared_ptr<double> J_ptr(0);
692 auto dJ_dX_vec = vectorDuplicate(direction_vec);
693 CHKERR VecZeroEntries(dJ_dX_vec);
694
695 auto fe_material =
696 boost::make_shared<VolumeElementForcesAndSourcesCore>(ep.mField);
697 CHKERR pushTopologicalMaterialOps(ep, fe_material, interior_integration_hook,
698 boundary_integration_hook, J_ptr,
699 dJ_dX_vec, eval_energy_model);
700 CHKERR DMoFEMLoopFiniteElements(ep.dmMaterial, ep.elementVolumeName,
701 fe_material);
702
703 CHKERR VecAssemblyBegin(dJ_dX_vec);
704 CHKERR VecAssemblyEnd(dJ_dX_vec);
705 CHKERR VecGhostUpdateBegin(dJ_dX_vec, ADD_VALUES, SCATTER_REVERSE);
706 CHKERR VecGhostUpdateEnd(dJ_dX_vec, ADD_VALUES, SCATTER_REVERSE);
707 CHKERR VecGhostUpdateBegin(dJ_dX_vec, INSERT_VALUES, SCATTER_FORWARD);
708 CHKERR VecGhostUpdateEnd(dJ_dX_vec, INSERT_VALUES, SCATTER_FORWARD);
709
710 auto delta_vec = vectorDuplicate(direction_vec);
711 CHKERR VecCopy(direction_vec, delta_vec);
712 CHKERR VecScale(delta_vec, epsilon);
713
714 auto a_vec = vectorDuplicate(delta_vec);
715 CHKERR VecCopy(row_sol, a_vec);
716 CHKERR VecAXPY(a_vec, 1, delta_vec);
717 auto b_vec = vectorDuplicate(delta_vec);
718 CHKERR VecCopy(row_sol, b_vec);
719 CHKERR VecAXPY(b_vec, -1, delta_vec);
720
721 auto fe_material_plus =
722 boost::make_shared<VolumeElementForcesAndSourcesCore>(ep.mField);
723 auto J_plus = boost::make_shared<double>(0);
724 CHKERR pushTopologicalMaterialOps(ep, fe_material_plus,
725 interior_integration_hook,
726 boundary_integration_hook, J_plus, nullptr,
727 eval_energy_model);
728 auto fe_material_minus =
729 boost::make_shared<VolumeElementForcesAndSourcesCore>(ep.mField);
730 auto J_minus = boost::make_shared<double>(0);
731 CHKERR pushTopologicalMaterialOps(
732 ep, fe_material_minus, interior_integration_hook,
733 boundary_integration_hook, J_minus, nullptr, eval_energy_model);
734
735 CHKERR setMaterialConfigurationImpl(ep, a_vec);
736 CHKERR DMoFEMLoopFiniteElements(ep.dmMaterial, ep.elementVolumeName,
737 fe_material_plus);
738 CHKERR DMoFEMLoopFiniteElements(ep.dmMaterial, ep.naturalBcElement,
739 fe_material_plus);
740 MOFEM_LOG("EP", Sev::inform)
741 << "Objective at a_vec: J_plus = " << std::setprecision(12) << *J_plus;
742
743 CHKERR setMaterialConfigurationImpl(ep, b_vec);
744 CHKERR DMoFEMLoopFiniteElements(ep.dmMaterial, ep.elementVolumeName,
745 fe_material_minus);
746 CHKERR DMoFEMLoopFiniteElements(ep.dmMaterial, ep.naturalBcElement,
747 fe_material_minus);
748 MOFEM_LOG("EP", Sev::inform)
749 << "Objective at b_vec: J_minus = " << std::setprecision(12) << *J_minus;
750
751 double dJ_da = (*J_plus - *J_minus) / (2 * epsilon);
752
753 double exact_dJ = 0;
754 CHKERR VecDot(dJ_dX_vec, direction_vec, &exact_dJ);
755 double fd_error = dJ_da - exact_dJ;
756
757 MOFEM_LOG("EP", Sev::inform)
758 << "Fd dJ/dX = " << dJ_da << ", exact dJ/dX = " << exact_dJ
759 << ", error = " << fd_error;
760
761 CHKERR setMaterialConfigurationImpl(ep, row_sol);
762 CHKERR setSpatialConfigurationImpl(ep, col_sol);
763
764 MoFEMFunctionReturn(0);
765}

◆ finiteDifference_dJdx_Test

MoFEMErrorCode finiteDifference_dJdx_Test ( TopologicalTAOCtxImpl *  ctx_impl_ptr,
double  epsilon,
ObjectiveModelType  eval_energy_model 
)
friend

Definition at line 515 of file EshelbianTopologicalDerivative.cpp.

516 {
517 MoFEMFunctionBegin;
518
519 MOFEM_LOG("EP", Sev::inform)
520 << "Starting finite difference dJ_dx gradient test with epsilon = "
521 << epsilon;
522
523 auto &ep = *ctx_impl_ptr->ep_ptr;
524
525 auto sol = createDMVector(ctx_impl_ptr->ep_ptr->dmElastic, RowColData::ROW);
526 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, sol, INSERT_VALUES,
527 SCATTER_FORWARD, RowColData::ROW);
528 CHKERR VecGhostUpdateBegin(sol, INSERT_VALUES, SCATTER_FORWARD);
529 CHKERR VecGhostUpdateEnd(sol, INSERT_VALUES, SCATTER_FORWARD);
530
531 const std::array<double, 2> piola_range{{-1, 1}};
532 const std::array<double, 2> bubble_range{{-1, 1}};
533 const std::array<double, 2> spatial_l2_disp_range{{-1, 1}};
534 const std::array<double, 2> rot_axis_range{{-1, 1}};
535 const std::array<double, 2> stretch_tensor_range{{-1, 1}};
536 const std::array<double, 2> hybrid_spatial_disp_range{{-1, 1}};
537
538 std::vector<OperatorsTester::RandomFieldData> random_fields{
539 {ep.piolaStress, piola_range},
540 {ep.bubbleField, bubble_range},
541 {ep.spatialL2Disp, spatial_l2_disp_range},
542 {ep.rotAxis, rot_axis_range},
543 {ep.stretchTensor, stretch_tensor_range},
544 {ep.hybridSpatialDisp, hybrid_spatial_disp_range}};
545
546 auto opt = ep.mField.getInterface<OperatorsTester>();
547 auto direction_vec = opt->setRandomFields(ep.dmElastic, random_fields,
548 nullptr, RowColData::ROW);
549
550 auto delta_vec = vectorDuplicate(direction_vec);
551 CHKERR VecCopy(direction_vec, delta_vec);
552 CHKERR VecScale(delta_vec, epsilon);
553
554 auto a_vec = vectorDuplicate(delta_vec);
555 CHKERR VecCopy(sol, a_vec);
556 CHKERR VecAXPY(a_vec, 1, delta_vec);
557 auto b_vec = vectorDuplicate(delta_vec);
558 CHKERR VecCopy(sol, b_vec);
559 CHKERR VecAXPY(b_vec, -1, delta_vec);
560
561 CHKERR setSpatialStateImpl(ep, a_vec);
562 double J_plus = 0;
563 CHKERR evaluateObjectiveImpl(ctx_impl_ptr, &J_plus, eval_energy_model);
564 MOFEM_LOG("EP", Sev::inform)
565 << "Objective at a_vec: J_plus = " << std::setprecision(12) << J_plus;
566
567 CHKERR setSpatialStateImpl(ep, b_vec);
568 double J_minus = 0;
569 CHKERR evaluateObjectiveImpl(ctx_impl_ptr, &J_minus, eval_energy_model);
570
571 MOFEM_LOG("EP", Sev::inform)
572 << "Objective at b_vec: J_minus = " << std::setprecision(12) << J_minus;
573 double dJ_da = (J_plus - J_minus) / (2 * epsilon);
574
575 double nrm_sol;
576 CHKERR VecNorm(sol, NORM_2, &nrm_sol);
577
578 CHKERR setSpatialStateImpl(ep, sol);
579
580 auto fe_spatial =
581 boost::make_shared<VolumeElementForcesAndSourcesCore>(ep.mField);
582 auto interior_integration_hook = ctx_impl_ptr->integrationAtInterior;
583 auto boundary_integration_hook = ctx_impl_ptr->integrationAtFace;
584
585 auto dJ_dx_vec =
586 pushTopologicalSpatialOps(ep, fe_spatial, interior_integration_hook,
587 boundary_integration_hook, eval_energy_model);
588 CHKERR VecZeroEntries(dJ_dx_vec);
589 CHKERR VecGhostUpdateBegin(dJ_dx_vec, INSERT_VALUES, SCATTER_FORWARD);
590 CHKERR VecGhostUpdateEnd(dJ_dx_vec, INSERT_VALUES, SCATTER_FORWARD);
591 CHKERR DMoFEMLoopFiniteElements(ep.dmElastic, ep.elementVolumeName,
592 fe_spatial);
593
594 CHKERR VecAssemblyBegin(dJ_dx_vec);
595 CHKERR VecAssemblyEnd(dJ_dx_vec);
596 CHKERR VecGhostUpdateBegin(dJ_dx_vec, ADD_VALUES, SCATTER_REVERSE);
597 CHKERR VecGhostUpdateEnd(dJ_dx_vec, ADD_VALUES, SCATTER_REVERSE);
598 CHKERR VecGhostUpdateBegin(dJ_dx_vec, INSERT_VALUES, SCATTER_FORWARD);
599 CHKERR VecGhostUpdateEnd(dJ_dx_vec, INSERT_VALUES, SCATTER_FORWARD);
600
601 double exact_dJ = 0;
602 CHKERR VecDot(dJ_dx_vec, direction_vec, &exact_dJ);
603 double error = dJ_da - exact_dJ;
604
605 MOFEM_LOG("EP", Sev::inform)
606 << "dJ/dx = " << dJ_da << ", exact dJ/dx = " << exact_dJ
607 << ", error = " << error << ", fraction " << dJ_da / exact_dJ
608 << " x norm: " << nrm_sol;
609
610
611
612
613 MoFEMFunctionReturn(0);
614}
EshelbianPlasticity::setSpatialStateImpl
static MoFEMErrorCode setSpatialStateImpl(EshelbianCore &ep, Vec sol)
Definition EshelbianTopologicalDerivative.cpp:505
EshelbianPlasticity::pushTopologicalSpatialOps
SmartPetscObj< Vec > pushTopologicalSpatialOps(EshelbianCore &ep, boost::shared_ptr< VolumeElementForcesAndSourcesCore > fe, ForcesAndSourcesCore::GaussHookFun interior_integration_hook, ForcesAndSourcesCore::GaussHookFun boundary_integration_hook, ObjectiveModelType eval_energy_model)
Definition EshelbianTopologicalDerivativePush.cpp:205
EshelbianPlasticity::TopologicalTAOCtxImpl::evaluateObjectiveImpl
friend MoFEMErrorCode evaluateObjectiveImpl(TopologicalTAOCtxImpl *ctx_impl_ptr, double *f, ObjectiveModelType eval_energy_model)
Definition EshelbianTopologicalDerivative.cpp:258

◆ finiteDifferenceGradientTest

MoFEMErrorCode finiteDifferenceGradientTest ( TopologicalTAOCtxImpl *  ctx_impl_ptr,
Vec  sol,
double *  f,
Vec  g,
double  epsilon,
ObjectiveModelType  eval_energy_model 
)
friend

Definition at line 317 of file EshelbianTopologicalDerivative.cpp.

320 {
321 MoFEMFunctionBegin;
322
323 MOFEM_LOG("EP", Sev::inform)
324 << "Starting finite difference dJ_dX gradient test with epsilon = "
325 << epsilon;
326
327 auto &ep = *ctx_impl_ptr->ep_ptr;
328 auto opt = ep.mField.getInterface<OperatorsTester>();
329
330 Range body_ents;
331 CHKERR ep.mField.get_moab().get_entities_by_dimension(0, SPACE_DIM,
332 body_ents);
333 auto body_skin = filter_true_skin(ep.mField, get_skin(ep.mField, body_ents));
334 Range body_skin_verts;
335 CHKERR ep.mField.get_moab().get_connectivity(body_skin, body_skin_verts,
336 true);
337
338 auto direction_vec = opt->setRandomFields(
339 ep.dmMaterial, {{ep.materialH1Positions, {-1, 1.}}},
340 boost::make_shared<Range>(body_skin_verts), RowColData::ROW);
341
342 auto delta_vec = vectorDuplicate(direction_vec);
343 CHKERR VecCopy(direction_vec, delta_vec);
344 CHKERR VecScale(delta_vec, epsilon);
345
346 auto a_vec = vectorDuplicate(delta_vec);
347 CHKERR VecCopy(sol, a_vec);
348 CHKERR VecAXPY(a_vec, 1., delta_vec);
349 auto b_vec = vectorDuplicate(delta_vec);
350 CHKERR VecCopy(sol, b_vec);
351 CHKERR VecAXPY(b_vec, -1., delta_vec);
352
353 CHKERR setMaterialConfigurationImpl(ep, a_vec);
354 CHKERR solveEquilibriumStateTS(ctx_impl_ptr);
355 double J_plus = 0;
356 CHKERR evaluateObjectiveImpl(ctx_impl_ptr, &J_plus, eval_energy_model);
357 MOFEM_LOG("EP", Sev::inform)
358 << "Objective at a_vec: J_plus = " << std::setprecision(12) << J_plus;
359
360 CHKERR setMaterialConfigurationImpl(ep, b_vec);
361 CHKERR solveEquilibriumStateTS(ctx_impl_ptr);
362 double J_minus = 0;
363 CHKERR evaluateObjectiveImpl(ctx_impl_ptr, &J_minus, eval_energy_model);
364
365 MOFEM_LOG("EP", Sev::inform)
366 << "Objective at b_vec: J_minus = " << std::setprecision(12) << J_minus;
367 double dJ_da = (J_plus - J_minus) / (2 * epsilon);
368
369 double exact_dJ = 0;
370 CHKERR VecDot(g, direction_vec, &exact_dJ);
371 double error = dJ_da - exact_dJ;
372
373 MOFEM_LOG("EP", Sev::inform)
374 << "J = " << *f << ", dJ/dX = " << dJ_da << ", exact dJ/dX = " << exact_dJ
375 << ", error = " << error << ", fraction " << dJ_da / exact_dJ;
376
377 auto *adj_problem_ptr = getProblemPtr(ep.dmMaterial);
378 auto &adj_dofs =
379 adj_problem_ptr->getNumeredRowDofsPtr()->get<PetscGlobalIdx_mi_tag>();
380
381 auto g_duplicate_vec = vectorDuplicate(g);
382 CHKERR VecZeroEntries(g_duplicate_vec);
383
384 SmartPetscObj<IS> is;
385 CHKERR ep.mField.getInterface<ISManager>()->isCreateProblemFieldAndRank(
386 adj_problem_ptr->getName(), RowColData::ROW, ep.materialH1Positions, 0, 3,
387 is, &body_skin_verts);
388 IS is_raw;
389 CHKERR ISAllGather(is, &is_raw);
390 is = SmartPetscObj<IS>(is_raw, false);
391 PetscInt nb_dofs;
392 CHKERR ISGetSize(is, &nb_dofs);
393 const PetscInt *is_ptr;
394 CHKERR ISGetIndices(is, &is_ptr);
395
396 constexpr double procent = 0; /* % */
397 const int nb_dofs_comp = ceil(procent * nb_dofs / 100.);
398
399 auto get_vec_value = [&](Vec vec, PetscInt idx) {
400 double *array;
401 CHKERR VecGetArray(vec, &array);
402 double &value = array[idx];
403 CHKERR VecRestoreArray(vec, &array);
404 return value;
405 };
406
407 auto set_vec_value = [&](Vec vec, PetscInt idx, double value) {
408 MoFEMFunctionBegin;
409 double *array;
410 CHKERR VecGetArray(vec, &array);
411 array[idx] = value;
412 CHKERR VecRestoreArray(vec, &array);
413 MoFEMFunctionReturn(0);
414 };
415
416 for (auto i = 0, j = 0; i != nb_dofs; ++i) {
417 MOFEM_LOG("EP", Sev::inform)
418 << "Testing DOF " << i << " out of " << nb_dofs;
419 auto dof_it = adj_dofs.find(is_ptr[i]);
420 if (dof_it != adj_dofs.end()) {
421 auto idx = (*dof_it)->getPetscLocalDofIdx();
422 auto exact_dJ = get_vec_value(g, idx);
423 constexpr double epsilon = 1e-8;
424 if (std::abs(exact_dJ) < epsilon) {
425 MOFEM_LOG("EP", Sev::inform)
426 << "Skipping DOF index " << idx
427 << " due to small gradient value: " << exact_dJ;
428 continue;
429 }
430 if (j >= nb_dofs_comp) {
431 MOFEM_LOG("EP", Sev::inform)
432 << "Stopping finite difference test after " << nb_dofs_comp
433 << " DOFs, out of total " << nb_dofs;
434 break;
435 }
436 ++j;
437
438 CHKERR VecCopy(sol, direction_vec);
439 CHKERR set_vec_value(direction_vec, idx,
440 get_vec_value(direction_vec, idx) + epsilon);
441 CHKERR setMaterialConfigurationImpl(ep, direction_vec);
442 CHKERR solveEquilibriumStateTS(ctx_impl_ptr);
443 double J_plus = 0;
444 CHKERR evaluateObjectiveImpl(ctx_impl_ptr, &J_plus, eval_energy_model);
445
446 CHKERR VecCopy(sol, direction_vec);
447 CHKERR set_vec_value(direction_vec, idx,
448 get_vec_value(direction_vec, idx) - epsilon);
449 CHKERR setMaterialConfigurationImpl(ep, direction_vec);
450 CHKERR solveEquilibriumStateTS(ctx_impl_ptr);
451 double J_minus = 0;
452 CHKERR evaluateObjectiveImpl(ctx_impl_ptr, &J_minus, eval_energy_model);
453
454 double dJ_da = (J_plus - J_minus) / (2 * epsilon);
455 CHKERR set_vec_value(g_duplicate_vec, idx, dJ_da);
456
457 double error = dJ_da - exact_dJ;
458 MOFEM_LOG("EP", Sev::inform)
459 << "DOF index: " << idx << ", J = " << *f << ", dJ/dX = " << dJ_da
460 << ", exact dJ/dX = " << exact_dJ << ", error = " << error
461 << ", fraction " << dJ_da / exact_dJ;
462 }
463 }
464
465
466 CHKERR ISRestoreIndices(is, &is_ptr);
467 CHKERR setMaterialConfigurationImpl(ep, sol);
468
469 double *g_array, *g_duplicate_array;
470 CHKERR VecGetArray(g, &g_array);
471 CHKERR VecGetArray(g_duplicate_vec, &g_duplicate_array);
472 auto &moab = ep.mField.get_moab();
473 Tag th_g, th_fd_g;
474 double def_val[] = {0, 0, 0};
475 CHKERR moab.tag_get_handle("G", 3, MB_TYPE_DOUBLE, th_g,
476 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
477 CHKERR moab.tag_get_handle("G_fd", 3, MB_TYPE_DOUBLE, th_fd_g,
478 MB_TAG_CREAT | MB_TAG_SPARSE, &def_val);
479 for (auto &dof : adj_dofs) {
480 if (!dof->getHasLocalIndex())
481 continue;
482 auto idx = dof->getPetscLocalDofIdx();
483 auto ent = dof->getEnt();
484 auto coeff = dof->getDofCoeffIdx();
485 if (coeff == 0) {
486 CHKERR moab.tag_set_data(th_g, &ent, 1, &g_array[idx]);
487 CHKERR moab.tag_set_data(th_fd_g, &ent, 1, &g_duplicate_array[idx]);
488 }
489 }
490 CHKERR VecRestoreArray(g, &g_array);
491 CHKERR VecRestoreArray(g_duplicate_vec, &g_duplicate_array);
492
493 EntityHandle root_mesh = 0;
494 std::vector<Tag> tags_list{th_g, th_fd_g};
495 CHKERR moab.write_file("gradient_comparison.h5m", "MOAB",
496 "PARALLEL=WRITE_PART", &root_mesh, 1,
497 &*tags_list.begin(), tags_list.size());
498
499 CHKERR moab.tag_delete(th_g);
500 CHKERR moab.tag_delete(th_fd_g);
501
502 MoFEMFunctionReturn(0);
503}
SPACE_DIM
constexpr int SPACE_DIM
Definition child_and_parent.cpp:17
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:21
EigenMatrix::Vec
const FTensor::Tensor2< T, Dim, Dim > Vec
Definition MatrixFunction.hpp:64
EshelbianPlasticity::filter_true_skin
static auto filter_true_skin(MoFEM::Interface &m_field, Range &&skin)
Definition EshelbianTopologicalDerivative.cpp:299
EshelbianPlasticity::get_skin
static auto get_skin(MoFEM::Interface &m_field, Range body_ents)
Definition EshelbianTopologicalDerivative.cpp:310
MoFEM::getProblemPtr
auto getProblemPtr(DM dm)
get problem pointer from DM
Definition DMMoFEM.hpp:1182
EshelbianPlasticity::TopologicalTAOCtxImpl::solveEquilibriumStateTS
friend MoFEMErrorCode solveEquilibriumStateTS(TopologicalTAOCtxImpl *ctx_impl_ptr)
Definition EshelbianTopologicalDerivative.cpp:116
MoFEM::ISManager
Section manager is used to create indexes and sections.
Definition ISManager.hpp:23
MoFEM::PetscGlobalIdx_mi_tag
Definition TagMultiIndices.hpp:38
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93

◆ solveAdjointSensitivityKSP

MoFEMErrorCode solveAdjointSensitivityKSP ( TopologicalTAOCtxImpl *  ctx_impl_ptr,
ObjectiveModelType  eval_energy_model 
)
friend

Definition at line 152 of file EshelbianTopologicalDerivative.cpp.

154 {
155 MoFEMFunctionBegin;
156 auto &ep = *ctx_impl_ptr->ep_ptr;
157 auto &m_field = ep.mField;
158 auto interior_integration_hook = ctx_impl_ptr->integrationAtInterior;
159 auto boundary_integration_hook = ctx_impl_ptr->integrationAtFace;
160
161 auto fe_spatial =
162 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
163 auto dJ_dx_vec =
164 pushTopologicalSpatialOps(ep, fe_spatial, interior_integration_hook,
165 boundary_integration_hook, eval_energy_model);
166 CHKERR VecZeroEntries(dJ_dx_vec);
167 CHKERR VecGhostUpdateBegin(dJ_dx_vec, INSERT_VALUES, SCATTER_FORWARD);
168 CHKERR VecGhostUpdateEnd(dJ_dx_vec, INSERT_VALUES, SCATTER_FORWARD);
169 CHKERR DMoFEMLoopFiniteElements(ep.dmElastic, ep.elementVolumeName,
170 fe_spatial);
171 CHKERR VecAssemblyBegin(dJ_dx_vec);
172 CHKERR VecAssemblyEnd(dJ_dx_vec);
173 CHKERR VecGhostUpdateBegin(dJ_dx_vec, ADD_VALUES, SCATTER_REVERSE);
174 CHKERR VecGhostUpdateEnd(dJ_dx_vec, ADD_VALUES, SCATTER_REVERSE);
175 CHKERR VecGhostUpdateBegin(dJ_dx_vec, INSERT_VALUES, SCATTER_FORWARD);
176 CHKERR VecGhostUpdateEnd(dJ_dx_vec, INSERT_VALUES, SCATTER_FORWARD);
177
178 double nrm_dJ_dx;
179 CHKERR VecNorm(dJ_dx_vec, NORM_2, &nrm_dJ_dx);
180 MOFEM_LOG("EP", Sev::inform)
181 << "Adjoint problem: dJ/dx vector norm: " << nrm_dJ_dx;
182
183 auto ksp = snesGetKSP(tsGetSNES(ctx_impl_ptr->timeSolver));
184 CHKERR KSPSolveTranspose(ksp, dJ_dx_vec, ctx_impl_ptr->adjointProblemVec);
185 CHKERR VecGhostUpdateBegin(ctx_impl_ptr->adjointProblemVec, INSERT_VALUES,
186 SCATTER_FORWARD);
187 CHKERR VecGhostUpdateEnd(ctx_impl_ptr->adjointProblemVec, INSERT_VALUES,
188 SCATTER_FORWARD);
189
190 double nrm_lambda;
191 CHKERR VecNorm(ctx_impl_ptr->adjointProblemVec, NORM_2, &nrm_lambda);
192 MOFEM_LOG("EP", Sev::inform)
193 << "Adjoint problem solved, lambda vector norm: " << nrm_lambda;
194
195 CHKERR ep.postProcessResults(0, "adjoint_solution.h5m", dJ_dx_vec,
196 ctx_impl_ptr->adjointProblemVec);
197
198 MoFEMFunctionReturn(0);
199}
MoFEM::snesGetKSP
auto snesGetKSP(SNES snes)
Definition PetscSmartObj.hpp:416
MoFEM::tsGetSNES
auto tsGetSNES(TS ts)
Definition PetscSmartObj.hpp:410

◆ solveEquilibriumStateTS

MoFEMErrorCode solveEquilibriumStateTS ( TopologicalTAOCtxImpl *  ctx_impl_ptr)
friend

Definition at line 116 of file EshelbianTopologicalDerivative.cpp.

117 {
118 MoFEMFunctionBegin;
119 auto &ep = *ctx_impl_ptr->ep_ptr;
120 auto ts = ctx_impl_ptr->timeSolver;
121
122 CHKERR VecGhostUpdateBegin(ctx_impl_ptr->primalProblemVec, INSERT_VALUES,
123 SCATTER_FORWARD);
124 CHKERR VecGhostUpdateEnd(ctx_impl_ptr->primalProblemVec, INSERT_VALUES,
125 SCATTER_FORWARD);
126 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, ctx_impl_ptr->primalProblemVec,
127 INSERT_VALUES, SCATTER_FORWARD,
128 RowColData::ROW);
129
130 CHKERR TSSetFromOptions(ts);
131 CHKERR TSSetStepNumber(ts, 0);
132 CHKERR TSSetTime(ts, 0);
133 CHKERR TSSetSolution(ts, ctx_impl_ptr->primalProblemVec);
134 double dt;
135 CHKERR TSGetTimeStep(ts, &dt);
136 MOFEM_LOG("EP", Sev::inform)
137 << "solveEquilibriumStateTS: Time step dt: " << dt;
138 CHKERR TSSolve(ts, PETSC_NULLPTR);
139 CHKERR TSSetTimeStep(ts, dt);
140
141 CHKERR VecGhostUpdateBegin(ctx_impl_ptr->primalProblemVec, INSERT_VALUES,
142 SCATTER_FORWARD);
143 CHKERR VecGhostUpdateEnd(ctx_impl_ptr->primalProblemVec, INSERT_VALUES,
144 SCATTER_FORWARD);
145 CHKERR DMoFEMMeshToLocalVector(ep.dmElastic, ctx_impl_ptr->primalProblemVec,
146 INSERT_VALUES, SCATTER_REVERSE,
147 RowColData::ROW);
148
149 MoFEMFunctionReturn(0);
150}
dt
double dt
Definition heat_method.cpp:27

◆ testTopologicalDerivative

MoFEMErrorCode testTopologicalDerivative ( TopologicalTAOCtx *  ctx_ptr,
Vec  sol,
PetscReal *  f,
Vec  g,
ObjectiveModelType  eval_energy_model 
)
friend

Definition at line 994 of file EshelbianTopologicalDerivative.cpp.

996 {
997
998 MoFEMFunctionBegin;
999 auto *ctx_impl_ptr = dynamic_cast<TopologicalTAOCtxImpl *>(ctx_ptr);
1000 if (!ctx_impl_ptr)
1001 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1002 "Invalid context pointer type");
1003 auto &ep = *ctx_impl_ptr->ep_ptr;
1004
1005 CHKERR VecGhostUpdateBegin(sol, INSERT_VALUES, SCATTER_FORWARD);
1006 CHKERR VecGhostUpdateEnd(sol, INSERT_VALUES, SCATTER_FORWARD);
1007 CHKERR DMoFEMMeshToLocalVector(ep.dmMaterial, sol, INSERT_VALUES,
1008 SCATTER_REVERSE, RowColData::ROW);
1009
1010 CHKERR VecZeroEntries(g);
1011 *f = 0;
1012
1013 CHKERR solveEquilibriumStateTS(ctx_impl_ptr);
1014 double norm2_x;
1015 CHKERR VecNorm(ctx_impl_ptr->primalProblemVec, NORM_2, &norm2_x);
1016 MOFEM_LOG("EP", Sev::inform)
1017 << "solveEquilibriumStateTS: Norm of displacement vector: "
1018 << norm2_x;
1019 if (norm2_x < 1e-12) {
1020 MoFEMFunctionReturnHot(0);
1021 }
1022
1023 double eps = 1e-6; // testing der
1024 CHKERR ::PetscOptionsGetReal(nullptr, nullptr, "-fd_epsilon", &eps, nullptr);
1025
1026 // derivative of J over X
1027 CHKERR finiteDifference_dJdX_Test(ctx_impl_ptr, eps, eval_energy_model);
1028 // derivative of J over x = { stress, strain, etc. }
1029 CHKERR finiteDifference_dJdx_Test(ctx_impl_ptr, eps, eval_energy_model);
1030 // derivative of J over X = { material distribution }
1031 CHKERR finiteDifference_dJd_adjoint_Test(ctx_impl_ptr, eps);
1032 // test dJ/dX = dJ/dX - dJ_adjoint/dX
1033
1034 // calculate lambda = K^ (-T) dJ/dx
1035 CHKERR solveAdjointSensitivityKSP(ctx_impl_ptr, eval_energy_model);
1036 // calculate dJ/dX = dJ/dX - dJ_adjoint/dX
1037 CHKERR evaluateGradientImpl(ctx_impl_ptr, f, g, eval_energy_model);
1038
1039 CHKERR ep.postProcessResults(0, "exact_gradient.h5m", PETSC_NULLPTR,
1040 PETSC_NULLPTR, g);
1041
1042 CHKERR finiteDifferenceGradientTest(ctx_impl_ptr, sol, f, g, eps,
1043 eval_energy_model);
1044
1045 double norm2_g; //< norm of the gradient vector
1046 CHKERR VecNorm(g, NORM_2, &norm2_g);
1047 MOFEM_LOG("EP", Sev::inform) << "Evaluated dJ_dX = " << norm2_g;
1048
1049
1050 MoFEMFunctionReturn(0);
1051}
eps
static const double eps
Definition check_base_functions_derivatives_on_tet.cpp:11
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
EshelbianPlasticity::TopologicalTAOCtxImpl::solveAdjointSensitivityKSP
friend MoFEMErrorCode solveAdjointSensitivityKSP(TopologicalTAOCtxImpl *ctx_impl_ptr, ObjectiveModelType eval_energy_model)
Definition EshelbianTopologicalDerivative.cpp:152
EshelbianPlasticity::TopologicalTAOCtxImpl::finiteDifference_dJdX_Test
friend MoFEMErrorCode finiteDifference_dJdX_Test(TopologicalTAOCtxImpl *ctx_impl_ptr, double epsilon, ObjectiveModelType eval_energy_model)
Definition EshelbianTopologicalDerivative.cpp:616
EshelbianPlasticity::TopologicalTAOCtxImpl::evaluateGradientImpl
friend MoFEMErrorCode evaluateGradientImpl(TopologicalTAOCtxImpl *ctx_impl_ptr, double *f, Vec g, ObjectiveModelType eval_energy_model)
Definition EshelbianTopologicalDerivative.cpp:201
EshelbianPlasticity::TopologicalTAOCtxImpl::finiteDifference_dJdx_Test
friend MoFEMErrorCode finiteDifference_dJdx_Test(TopologicalTAOCtxImpl *ctx_impl_ptr, double epsilon, ObjectiveModelType eval_energy_model)
Definition EshelbianTopologicalDerivative.cpp:515
EshelbianPlasticity::TopologicalTAOCtxImpl::finiteDifferenceGradientTest
friend MoFEMErrorCode finiteDifferenceGradientTest(TopologicalTAOCtxImpl *ctx_impl_ptr, Vec sol, double *f, Vec g, double epsilon, ObjectiveModelType eval_energy_model)
Definition EshelbianTopologicalDerivative.cpp:317
EshelbianPlasticity::TopologicalTAOCtxImpl::finiteDifference_dJd_adjoint_Test
friend MoFEMErrorCode finiteDifference_dJd_adjoint_Test(TopologicalTAOCtxImpl *ctx_impl_ptr, double epsilon)
Definition EshelbianTopologicalDerivative.cpp:768

◆ topologicalEvaluateObjectiveAndGradient [1/2]

MoFEMErrorCode topologicalEvaluateObjectiveAndGradient ( Tao  tao,
Vec  sol,
PetscReal *  f,
Vec  g,
void *  ctx 
)
friend

Definition at line 1071 of file EshelbianTopologicalDerivative.cpp.

1073 {
1074
1075 MoFEMFunctionBegin;
1076 // auto *ctx_impl_ptr = static_cast<TopologicalTAOCtxImpl *>(ctx);
1077 // if (!ctx_impl_ptr)
1078 // SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1079 // "Invalid context pointer type");
1080 // auto &ep = *ctx_impl_ptr->ep_ptr;
1081 // auto &m_field = ep.mField;
1082 // auto interior_integration_hook = ctx_impl_ptr->integrationAtInterior;
1083 // auto boundary_integration_hook = ctx_impl_ptr->integrationAtFace;
1084 // auto alpha = ep.alphaW;
1085 // auto rho = ep.alphaRho;
1086 // auto alpha_omega = ep.alphaOmega;
1087
1088 MoFEMFunctionReturn(0);
1089}

◆ topologicalEvaluateObjectiveAndGradient [2/2]

MoFEMErrorCode topologicalEvaluateObjectiveAndGradient ( Tao  tao,
Vec  sol,
PetscReal *  f,
Vec  g,
void *  ctx 
)
friend

Definition at line 1071 of file EshelbianTopologicalDerivative.cpp.

1073 {
1074
1075 MoFEMFunctionBegin;
1076 // auto *ctx_impl_ptr = static_cast<TopologicalTAOCtxImpl *>(ctx);
1077 // if (!ctx_impl_ptr)
1078 // SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1079 // "Invalid context pointer type");
1080 // auto &ep = *ctx_impl_ptr->ep_ptr;
1081 // auto &m_field = ep.mField;
1082 // auto interior_integration_hook = ctx_impl_ptr->integrationAtInterior;
1083 // auto boundary_integration_hook = ctx_impl_ptr->integrationAtFace;
1084 // auto alpha = ep.alphaW;
1085 // auto rho = ep.alphaRho;
1086 // auto alpha_omega = ep.alphaOmega;
1087
1088 MoFEMFunctionReturn(0);
1089}

Member Data Documentation

◆ adjointProblemVec

SmartPetscObj<Vec> EshelbianPlasticity::TopologicalTAOCtxImpl::adjointProblemVec
private

Definition at line 98 of file EshelbianTopologicalDerivative.cpp.

◆ ep_ptr

EshelbianCore* EshelbianPlasticity::TopologicalTAOCtxImpl::ep_ptr
private

Definition at line 46 of file EshelbianTopologicalDerivative.cpp.

◆ integrationAtFace

ForcesAndSourcesCore::GaussHookFun EshelbianPlasticity::TopologicalTAOCtxImpl::integrationAtFace
private

Definition at line 48 of file EshelbianTopologicalDerivative.cpp.

◆ integrationAtInterior

ForcesAndSourcesCore::GaussHookFun EshelbianPlasticity::TopologicalTAOCtxImpl::integrationAtInterior
private

Definition at line 47 of file EshelbianTopologicalDerivative.cpp.

◆ primalProblemVec

SmartPetscObj<Vec> EshelbianPlasticity::TopologicalTAOCtxImpl::primalProblemVec
private

Definition at line 97 of file EshelbianTopologicalDerivative.cpp.

◆ timeSolver

SmartPetscObj<TS> EshelbianPlasticity::TopologicalTAOCtxImpl::timeSolver
private

Definition at line 49 of file EshelbianTopologicalDerivative.cpp.

The documentation for this struct was generated from the following file:
Generated by Doxygen 1.9.8 and hosted at University of Glasgow