v0.16.0
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 evaluatePrimalProblemTopologicalImpl (TopologicalTAOCtxImpl *ctx_impl_ptr)
 
MoFEMErrorCode evaluateAdjointProblemImpl (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 31 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 32 of file EshelbianTopologicalDerivative.cpp.

37 : ep_ptr(ep), integrationAtInterior(set_integration_at_interior),
38 integrationAtFace(set_integration_at_face), timeSolver(time_solver) {
39
40 primalProblemVec = createDMVector(ep->dmElastic, RowColData::COL);
41 adjointProblemVec = createDMVector(ep->dmElastic, RowColData::COL);
42 }
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:270

Friends And Related Symbol Documentation

◆ evaluateAdjointProblemImpl

MoFEMErrorCode evaluateAdjointProblemImpl ( TopologicalTAOCtxImpl *  ctx_impl_ptr,
ObjectiveModelType  eval_energy_model 
)
friend

Definition at line 151 of file EshelbianTopologicalDerivative.cpp.

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

◆ evaluateGradientImpl

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

Definition at line 199 of file EshelbianTopologicalDerivative.cpp.

201 {
202 MoFEMFunctionBegin;
203 auto &ep = *ctx_impl_ptr->ep_ptr;
204 auto &m_field = ep.mField;
205 auto interior_integration_hook = ctx_impl_ptr->integrationAtInterior;
206 auto boundary_integration_hook = ctx_impl_ptr->integrationAtFace;
207
208 auto fe_material =
209 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
210 auto fe_interior_adjoint =
211 boost::make_shared<VolumeElementForcesAndSourcesCore>(m_field);
212 auto fe_natural_adjoint =
213 boost::make_shared<FaceElementForcesAndSourcesCore>(m_field);
214
215 boost::shared_ptr<double> J_ptr(
216 f, [](double *) {}); // Custom deleter to avoid deallocation
217 auto dJ_dX_vec = SmartPetscObj<Vec>(g, true);
218 CHKERR VecZeroEntries(dJ_dX_vec);
219 CHKERR VecGhostUpdateBegin(dJ_dX_vec, INSERT_VALUES, SCATTER_FORWARD);
220 CHKERR VecGhostUpdateEnd(dJ_dX_vec, INSERT_VALUES, SCATTER_FORWARD);
221
222 CHKERR pushTopologicalMaterialOps(ep, fe_material, interior_integration_hook,
223 boundary_integration_hook, J_ptr,
224 dJ_dX_vec, eval_energy_model);
225
226 CHKERR VecScale(ctx_impl_ptr->adjointProblemVec, -1.0);
227 auto alpha = ep.alphaW;
228 auto rho = ep.alphaRho;
229 auto alpha_omega = ep.alphaOmega;
230 CHKERR pushTopologicalInteriorOps_dJ_adjoint_gradient(
231 ep, fe_interior_adjoint, interior_integration_hook,
232 boundary_integration_hook, ctx_impl_ptr->adjointProblemVec, dJ_dX_vec,
233 alpha, rho, alpha_omega, nullptr);
234 CHKERR pushTopologicalBoundaryOps_dJ_adjoint_gradient(
235 ep, fe_natural_adjoint, interior_integration_hook,
236 boundary_integration_hook, ctx_impl_ptr->adjointProblemVec, dJ_dX_vec,
237 nullptr);
238
239 CHKERR DMoFEMLoopFiniteElements(ep.dmMaterial, ep.elementVolumeName,
240 fe_material);
241 CHKERR DMoFEMLoopFiniteElements(ep.dmMaterial, ep.elementVolumeName,
242 fe_interior_adjoint);
243 CHKERR DMoFEMLoopFiniteElements(ep.dmMaterial, ep.naturalBcElement,
244 fe_natural_adjoint);
245
246 CHKERR VecAssemblyBegin(dJ_dX_vec);
247 CHKERR VecAssemblyEnd(dJ_dX_vec);
248 CHKERR VecGhostUpdateBegin(dJ_dX_vec, ADD_VALUES, SCATTER_REVERSE);
249 CHKERR VecGhostUpdateEnd(dJ_dX_vec, ADD_VALUES, SCATTER_REVERSE);
250 CHKERR VecGhostUpdateBegin(dJ_dX_vec, INSERT_VALUES, SCATTER_FORWARD);
251 CHKERR VecGhostUpdateEnd(dJ_dX_vec, INSERT_VALUES, SCATTER_FORWARD);
252
253 MoFEMFunctionReturn(0);
254}
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:463
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:259
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_omega, boost::shared_ptr< double > J_ptr, SmartPetscObj< Vec > topo_vec=SmartPetscObj< Vec >())
Definition EshelbianTopologicalDerivativePush.cpp:314
g
constexpr double g
Definition shallow_wave.cpp:63
MoFEM::SmartPetscObj
intrusive_ptr for managing petsc objects
Definition PetscSmartObj.hpp:85
rho
double rho
Definition plastic.cpp:145

◆ evaluateObjectiveImpl

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

Definition at line 256 of file EshelbianTopologicalDerivative.cpp.

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

◆ evaluatePrimalProblemTopologicalImpl

MoFEMErrorCode evaluatePrimalProblemTopologicalImpl ( TopologicalTAOCtxImpl *  ctx_impl_ptr)
friend

Definition at line 116 of file EshelbianTopologicalDerivative.cpp.

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

◆ finiteDifference_dJd_adjoint_Test

MoFEMErrorCode finiteDifference_dJd_adjoint_Test ( TopologicalTAOCtxImpl *  ctx_impl_ptr,
double  epsilon 
)
friend

Definition at line 766 of file EshelbianTopologicalDerivative.cpp.

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

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

◆ finiteDifference_dJdx_Test

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

Definition at line 513 of file EshelbianTopologicalDerivative.cpp.

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

◆ finiteDifferenceGradientTest

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

Definition at line 315 of file EshelbianTopologicalDerivative.cpp.

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

◆ testTopologicalDerivative

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

Definition at line 992 of file EshelbianTopologicalDerivative.cpp.

994 {
995
996 MoFEMFunctionBegin;
997 auto *ctx_impl_ptr = dynamic_cast<TopologicalTAOCtxImpl *>(ctx_ptr);
998 if (!ctx_impl_ptr)
999 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1000 "Invalid context pointer type");
1001 auto &ep = *ctx_impl_ptr->ep_ptr;
1002
1003 CHKERR VecGhostUpdateBegin(sol, INSERT_VALUES, SCATTER_FORWARD);
1004 CHKERR VecGhostUpdateEnd(sol, INSERT_VALUES, SCATTER_FORWARD);
1005 CHKERR DMoFEMMeshToLocalVector(ep.dmMaterial, sol, INSERT_VALUES,
1006 SCATTER_REVERSE, RowColData::ROW);
1007
1008 CHKERR VecZeroEntries(g);
1009 *f = 0;
1010
1011 CHKERR evaluatePrimalProblemTopologicalImpl(ctx_impl_ptr);
1012 double norm2_x;
1013 CHKERR VecNorm(ctx_impl_ptr->primalProblemVec, NORM_2, &norm2_x);
1014 MOFEM_LOG("EP", Sev::inform)
1015 << "evaluatePrimalProblemTopologicalImpl: Norm of displacement vector: "
1016 << norm2_x;
1017 if (norm2_x < 1e-12) {
1018 MoFEMFunctionReturnHot(0);
1019 }
1020
1021 double eps = 1e-6; // testing der
1022 CHKERR ::PetscOptionsGetReal(nullptr, nullptr, "-fd_epsilon", &eps, nullptr);
1023
1024 // derivative of J over X
1025 CHKERR finiteDifference_dJdX_Test(ctx_impl_ptr, eps, eval_energy_model);
1026 // derivative of J over x = { stress, strain, etc. }
1027 CHKERR finiteDifference_dJdx_Test(ctx_impl_ptr, eps, eval_energy_model);
1028 // derivative of J over X = { material distribution }
1029 CHKERR finiteDifference_dJd_adjoint_Test(ctx_impl_ptr, eps);
1030 // test dJ/dX = dJ/dX - dJ_adjoint/dX
1031
1032 // calculate lambda = K^ (-T) dJ/dx
1033 CHKERR evaluateAdjointProblemImpl(ctx_impl_ptr, eval_energy_model);
1034 // calculate dJ/dX = dJ/dX - dJ_adjoint/dX
1035 CHKERR evaluateGradientImpl(ctx_impl_ptr, f, g, eval_energy_model);
1036
1037 CHKERR ep.postProcessResults(0, "exact_gradient.h5m", PETSC_NULLPTR,
1038 PETSC_NULLPTR, g);
1039
1040 CHKERR finiteDifferenceGradientTest(ctx_impl_ptr, sol, f, g, eps,
1041 eval_energy_model);
1042
1043 double norm2_g; //< norm of the gradient vector
1044 CHKERR VecNorm(g, NORM_2, &norm2_g);
1045 MOFEM_LOG("EP", Sev::inform) << "Evaluated dJ_dX = " << norm2_g;
1046
1047
1048 MoFEMFunctionReturn(0);
1049}
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::finiteDifference_dJdX_Test
friend MoFEMErrorCode finiteDifference_dJdX_Test(TopologicalTAOCtxImpl *ctx_impl_ptr, double epsilon, ObjectiveModelType eval_energy_model)
Definition EshelbianTopologicalDerivative.cpp:614
EshelbianPlasticity::TopologicalTAOCtxImpl::evaluateGradientImpl
friend MoFEMErrorCode evaluateGradientImpl(TopologicalTAOCtxImpl *ctx_impl_ptr, double *f, Vec g, ObjectiveModelType eval_energy_model)
Definition EshelbianTopologicalDerivative.cpp:199
EshelbianPlasticity::TopologicalTAOCtxImpl::evaluateAdjointProblemImpl
friend MoFEMErrorCode evaluateAdjointProblemImpl(TopologicalTAOCtxImpl *ctx_impl_ptr, ObjectiveModelType eval_energy_model)
Definition EshelbianTopologicalDerivative.cpp:151
EshelbianPlasticity::TopologicalTAOCtxImpl::finiteDifference_dJdx_Test
friend MoFEMErrorCode finiteDifference_dJdx_Test(TopologicalTAOCtxImpl *ctx_impl_ptr, double epsilon, ObjectiveModelType eval_energy_model)
Definition EshelbianTopologicalDerivative.cpp:513
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:315
EshelbianPlasticity::TopologicalTAOCtxImpl::finiteDifference_dJd_adjoint_Test
friend MoFEMErrorCode finiteDifference_dJd_adjoint_Test(TopologicalTAOCtxImpl *ctx_impl_ptr, double epsilon)
Definition EshelbianTopologicalDerivative.cpp:766

◆ topologicalEvaluateObjectiveAndGradient [1/2]

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

Definition at line 1069 of file EshelbianTopologicalDerivative.cpp.

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

◆ topologicalEvaluateObjectiveAndGradient [2/2]

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

Definition at line 1069 of file EshelbianTopologicalDerivative.cpp.

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

Member Data Documentation

◆ adjointProblemVec

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

Definition at line 97 of file EshelbianTopologicalDerivative.cpp.

◆ ep_ptr

EshelbianCore* EshelbianPlasticity::TopologicalTAOCtxImpl::ep_ptr
private

Definition at line 45 of file EshelbianTopologicalDerivative.cpp.

◆ integrationAtFace

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

Definition at line 47 of file EshelbianTopologicalDerivative.cpp.

◆ integrationAtInterior

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

Definition at line 46 of file EshelbianTopologicalDerivative.cpp.

◆ primalProblemVec

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

Definition at line 96 of file EshelbianTopologicalDerivative.cpp.

◆ timeSolver

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

Definition at line 48 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