Skip to first unread message
Assigned to lik...@wp.pl by me
Karol Lewandowski
Dec 6, 2021, 7:32:13 AM12/6/21
to MoFEM Q&A
Hello,
I am trying to create a subproblem for contact, I generate it in the following way:
MoFEMErrorCode set_contact_dm(MoFEM::Interface &m_field, SmartPetscObj<DM> &dm,
SmartPetscObj<DM> &dm_contact) {
MoFEMFunctionBegin;
auto simple = m_field.getInterface<Simple>();
dm_contact = createSmartDM(m_field.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(dm_contact, dm, "CONTACT_PROBLEM");
CHKERR DMMoFEMSetDestroyProblem(dm_contact, PETSC_TRUE);
CHKERR DMMoFEMAddSubFieldRow(dm_contact, "U");
if (m_field.check_field("TAU"))
CHKERR DMMoFEMAddSubFieldRow(dm_contact, "TAU");
if (m_field.check_field("EP"))
CHKERR DMMoFEMAddSubFieldRow(dm_contact, "EP");
CHKERR DMMoFEMAddElement(dm_contact, simple->getDomainFEName().c_str());
CHKERR DMMoFEMAddElement(dm_contact, simple->getBoundaryFEName().c_str());
// CHKERR DMMoFEMAddElement(dm_contact, simple->getSkeletonFEName().c_str());
CHKERR DMMoFEMSetSquareProblem(dm_contact, PETSC_TRUE);
CHKERR DMMoFEMSetIsPartitioned(dm_contact, PETSC_TRUE);
CHKERR DMSetUp(dm_contact);
if (1) {
SmartPetscObj<Mat> m;
CHKERR m_field.getInterface<MatrixManager>()
->createMPIAIJWithArrays<PetscGlobalIdx_mi_tag>("CONTACT_PROBLEM", m);
// MatView(m, PETSC_VIEWER_STDOUT_WORLD);
// MatView(m, PETSC_VIEWER_DRAW_WORLD);
// MatView(m, PETSC_VIEWER_MATLAB_WORLD);
// std::string wait;
// std::cin >> wait;
CHKERR m_field.getInterface<MatrixManager>()
->checkMPIAIJWithArraysMatrixFillIn<PetscGlobalIdx_mi_tag>(
"CONTACT_PROBLEM", 1, 1, 1);
}
MoFEMFunctionReturn(0);
}
As you can see at the end I run "checkMPIAIJWithArraysMatrixFillIn" and get the following warnings:
[0] <warning> Warning: Number of Dofs in Col diffrent than number of dofs for given entity order 24 28[0] <warning> Warning: Number of Dofs in Row diffrent than number of dofs for given entity order 4 28
Any ideas of what am I doing wrong? It seemed straightforwards
Karol Lewandowski
Dec 6, 2021, 7:34:06 AM12/6/21
to MoFEM Q&A
And later I have a problem when I try to assemble some fields. The error is:
```
0]PETSC ERROR: --------------------- Error Message --------------------------------------------------------------
[0]PETSC ERROR: Argument out of range
[0]PETSC ERROR: Inserting a new nonzero at global row/column (942, 944) into matrix
[0]PETSC ERROR: See http://www.mcs.anl.gov/petsc/documentation/faq.html for trouble shooting.
[0]PETSC ERROR: Petsc Release Version 3.11.4, Sep, 28, 2019
```
Message has been deleted
Lukasz Kaczmraczyk
Dec 6, 2021, 9:41:44 AM12/6/21
to MoFEM Q&A
Hello,
Thank you for your question!
I do not see anything wrong. The warning suggests that on a given row/col, i.e. entity associated with it is more DOFs, you had removed DOFs from a problem in DM. So those entries are not available in sub dm and associated problem "CONTACT_PROBLEM".
That is not an issue by itself but an unusual case; that is why it is a warning since is you typically test the consistency of problem structures before you strat remove DOFs.
In terms of assembly, all is OK; if that is an issue, you will get an error.
I think an error later is rather caused by something else, i.e. wrong indices. Since checkMPIAIJWithArraysMatrixFillIn makes crude test, i.e. it trying to assemble every possible entry given on element. So if you have an error, you try to assemble something that is not given by data.getIndices().
I hope that this helps.
Regards,
Lukasz
Karol Lewandowski
Dec 6, 2021, 7:26:55 PM12/6/21
to MoFEM Q&A
So the thing is that everything is ok with Indices until I add the contact on top, somehow I am not passing the index data correctly, even though I created AO.
Maybe you can spot some obvious mistake:
// Setup fieldsplit (block) solver - optional: yes/no
if (is_pcfs == PETSC_TRUE && !data.is_fieldsplit_bc_only) {
CHKERR set_global_mat_and_vec();
PetscBool is_l0_field_split; // level0 fieldsplit
PetscBool is_l1_field_split = PETSC_TRUE; // level1 fieldsplit
PetscObjectTypeCompare((PetscObject)pC, PCFIELDSPLIT, &is_l0_field_split);
// optional fieldsplit on top
auto set_fieldsplit_contact = [&]() {
MoFEMFunctionBeginHot;
const MoFEM::Problem *fs_sig_ptr;
CHKERR DMMoFEMGetProblemPtr(dmContact, &fs_sig_ptr);
if (auto sig_data = fs_sig_ptr->subProblemData) {
is_map["E_IS_SIG"] = sig_data->rowIs;
sigma_ao = sig_data->getSmartRowMap();
CHKERR PCFieldSplitSetIS(pC, NULL, is_map["SIGMA"]);
CHKERR PCFieldSplitSetIS(pC, NULL, is_map["E_IS_SIG"]);
CHKERR PCSetUp(pC);
PetscInt n;
KSP *ct_ksp;
CHKERR PCFieldSplitGetSubKSP(pC, &n, &ct_ksp);
PC l1_pc;
CHKERR KSPGetPC(ct_ksp[1], &l1_pc);
PetscObjectTypeCompare((PetscObject)l1_pc, PCFIELDSPLIT,
&is_l1_field_split);
// important!
// in case of contact analysis we swap the first preconditioner
pC = l1_pc;
CHKERR PetscFree(ct_ksp);
}
MoFEMFunctionReturnHot(0);
};
if (data.with_contact && is_l0_field_split)
CHKERR set_fieldsplit_contact();
PetscObjectTypeCompare((PetscObject)pC, PCFIELDSPLIT, &is_l1_field_split);
if (data.with_plasticity && is_l0_field_split && is_l1_field_split) {
const MoFEM::Problem *schur_ep_ptr;
CHKERR DMMoFEMGetProblemPtr(dmEP, &schur_ep_ptr);
if (auto ep_data = schur_ep_ptr->subProblemData) {
is_map["E_IS_EP"] = ep_data->rowIs;
if (sigma_ao)
CHKERR AOApplicationToPetscIS(sigma_ao, is_map["EP"]);
CHKERR PCFieldSplitSetIS(pC, NULL, is_map["EP"]);
CHKERR PCFieldSplitSetIS(pC, NULL, is_map["E_IS_EP"]);
PCCompositeType pc_type;
CHKERR PCFieldSplitGetType(pC, &pc_type);
if (pc_type == PC_COMPOSITE_SCHUR) {
CHKERR DMCreateMatrix_MoFEM(dmEP, commonDataPtr->SEpEp);
commonDataPtr->aoSEpEp = ep_data->getSmartRowMap();
commonDataPtr->blockMatContainerPtr =
boost::make_shared<BlockMatContainer>();
MoFEM::block_mat_container_ptr =
commonDataPtr->blockMatContainerPtr;
CHKERR PCFieldSplitSetSchurPre(pC, PC_FIELDSPLIT_SCHUR_PRE_USER,
commonDataPtr->SEpEp);
CHKERR PCSetUp(pC);
}
PetscInt n;
KSP *tt_ksp;
CHKERR PCFieldSplitGetSubKSP(pC, &n, &tt_ksp);
PC tau_pc;
CHKERR KSPGetPC(tt_ksp[1], &tau_pc);
PetscBool is_tau_field_split;
PetscObjectTypeCompare((PetscObject)tau_pc, PCFIELDSPLIT,
&is_tau_field_split);
if (is_tau_field_split) {
const MoFEM::Problem *schur_tau_ptr;
CHKERR DMMoFEMGetProblemPtr(dmTAU, &schur_tau_ptr);
if (auto tau_data = schur_tau_ptr->subProblemData) {
// CHKERR tau_data->getRowIs(is_map["E_IS_TAU"]);
is_map["E_IS_TAU"] = tau_data->rowIs;
AO ep_ao;
CHKERR ep_data->getRowMap(&ep_ao);
if (sigma_ao)
CHKERR AOApplicationToPetscIS(sigma_ao, is_map["TAU"]);
CHKERR AOApplicationToPetscIS(ep_ao, is_map["TAU"]);
CHKERR PCFieldSplitSetIS(tau_pc, NULL, is_map["TAU"]);
CHKERR PCFieldSplitSetIS(tau_pc, NULL, is_map["E_IS_TAU"]);
CHKERR DMCreateMatrix_MoFEM(dmTAU, commonDataPtr->STauTau);
commonDataPtr->aoSTauTau = tau_data->getSmartRowMap();
CHKERR PCFieldSplitSetSchurPre(
tau_pc, PC_FIELDSPLIT_SCHUR_PRE_USER, commonDataPtr->STauTau);
CHKERR PCSetUp(tau_pc);
PetscInt bc_n;
KSP *bc_ksp;
CHKERR PCFieldSplitGetSubKSP(tau_pc, &bc_n, &bc_ksp);
PC bc_pc;
CHKERR KSPGetPC(bc_ksp[1], &bc_pc);
CHKERR set_fieldsplit_on_bc(bc_pc, schur_tau_ptr->getName());
CHKERR PetscFree(bc_ksp);
}
}
CHKERR PetscFree(tt_ksp);
// for (auto &m : is_map)
// CHKERR ISDestroy(&m.second);
}
}
}
Lukasz Kaczmraczyk
Dec 7, 2021, 7:20:56 AM12/7/21
to MoFEM Q&A
Karol,
I think, that problem is create at places like this,
if (sigma_ao)
CHKERR AOApplicationToPetscIS(sigma_ao, is_map["TAU"]);
Since sigma_ao, you can not apply AO from one level up first. However I do not fully comprehend code at that stage, all implications of it. You have to apply cascade of AO. I do that probably in Elshelbian plasticity.
I will suggest to do, to dissect problem. What will happen when:
1. you apply FS only to SIGMA, and rest solve with direct solver
2. you apply FS to SIGMA then to EP, and solve rest with direct solver
3. if you do FS SIGMA->EP->TAU.
I hope that helps.
Kind regards,
Lukasz
Karol Lewandowski
Dec 7, 2021, 7:34:29 AM12/7/21
to MoFEM Q&A
"Since sigma_ao, you can not apply AO from one level up first. However I do not fully comprehend code at that stage, all implications of it. You have to apply cascade of AO. I do that probably in Elshelbian plasticity."
What do you mean?
I apply consecutive AOs there, first SIGMA (it is AO smart pointer) and then EP, otherwise, it will not work at all. PC setup would fail.
Lukasz Kaczmraczyk
Dec 7, 2021, 8:19:03 AM12/7/21
to MoFEM Q&A
Hello,
In fact, when I look at this
if (sigma_ao)
CHKERR AOApplicationToPetscIS(sigma_ao, is_map["TAU"]);
CHKERR AOApplicationToPetscIS(ep_ao, is_map["TAU"]);
that looks OK.
To find problem, we have to work step by step. What happens when you apply one FS? Then nested, and finally double nested.
Regards,
Lukasz
Karol Lewandowski
Dec 7, 2021, 10:04:27 AM12/7/21
to MoFEM Q&A
When I apply a single FS contact (additive or multiplicative), it works. The second level with multiplicative (EP) - it works as well,
Second level with Schur - problems with assembly.
Lukasz Kaczmraczyk
Dec 7, 2021, 10:15:40 AM12/7/21
to MoFEM Q&A
Ok. That suggest that we have problems with AO-s.
Lukasz Kaczmraczyk
Dec 7, 2021, 10:17:33 AM12/7/21
to MoFEM Q&A
We missing how you construct sub dm for tau? Can you show this?
Karol Lewandowski
Dec 7, 2021, 10:29:14 AM12/7/21
to MoFEM Q&A
I think sub dm for tau is not the issue here, already for EP I am getting assembly errors. Here they are:
in the main, this is how I set them:
MoFEMErrorCode set_plastic_ep_dm(MoFEM::Interface &m_field,
SmartPetscObj<DM> &dm,
SmartPetscObj<DM> &dm_plastic) {
MoFEMFunctionBegin;
auto simple = m_field.getInterface<Simple>();
dm_plastic = createSmartDM(m_field.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(dm_plastic, dm, "PLASTIC_EP_PROBLEM");
CHKERR DMMoFEMSetDestroyProblem(dm_plastic, PETSC_TRUE);
CHKERR DMMoFEMAddSubFieldRow(dm_plastic, "U");
CHKERR DMMoFEMAddSubFieldRow(dm_plastic, "TAU");
// if (m_field.check_field("SIGMA"))
// CHKERR DMMoFEMAddSubFieldRow(dm_plastic, "SIGMA");
CHKERR DMMoFEMAddElement(dm_plastic, simple->getDomainFEName().c_str());
CHKERR DMMoFEMAddElement(dm_plastic, simple->getBoundaryFEName().c_str());
CHKERR DMMoFEMSetSquareProblem(dm_plastic, PETSC_TRUE);
CHKERR DMMoFEMSetIsPartitioned(dm_plastic, PETSC_TRUE);
CHKERR DMSetUp(dm_plastic);
MoFEMFunctionReturn(0);
}
MoFEMErrorCode set_plastic_tau_dm(MoFEM::Interface &m_field,
SmartPetscObj<DM> &dm,
SmartPetscObj<DM> &dm_plastic_tau) {
MoFEMFunctionBegin;
auto simple = m_field.getInterface<Simple>();
dm_plastic_tau = createSmartDM(m_field.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(dm_plastic_tau, dm, "PLASTIC_TAU_PROBLEM");
CHKERR DMMoFEMSetDestroyProblem(dm_plastic_tau, PETSC_TRUE);
CHKERR DMMoFEMAddSubFieldRow(dm_plastic_tau, "U");
// if (m_field.check_field("SIGMA"))
// CHKERR DMMoFEMAddSubFieldRow(dm_plastic_tau, "SIGMA");
CHKERR DMMoFEMAddElement(dm_plastic_tau, simple->getDomainFEName().c_str());
CHKERR DMMoFEMAddElement(dm_plastic_tau, simple->getBoundaryFEName().c_str());
CHKERR DMMoFEMSetSquareProblem(dm_plastic_tau, PETSC_TRUE);
CHKERR DMMoFEMSetIsPartitioned(dm_plastic_tau, PETSC_TRUE);
CHKERR DMSetUp(dm_plastic_tau);
MoFEMFunctionReturn(0);
}
MoFEMErrorCode set_contact_dm(MoFEM::Interface &m_field, SmartPetscObj<DM> &dm,
SmartPetscObj<DM> &dm_contact) {
MoFEMFunctionBegin;
auto simple = m_field.getInterface<Simple>();
dm_contact = createSmartDM(m_field.get_comm(), "DMMOFEM");
CHKERR DMMoFEMCreateSubDM(dm_contact, dm, "CONTACT_PROBLEM");
CHKERR DMMoFEMSetDestroyProblem(dm_contact, PETSC_TRUE);
CHKERR DMMoFEMAddSubFieldRow(dm_contact, "U");
if (m_field.check_field("TAU"))
CHKERR DMMoFEMAddSubFieldRow(dm_contact, "TAU");
if (m_field.check_field("EP"))
CHKERR DMMoFEMAddSubFieldRow(dm_contact, "EP");
CHKERR DMMoFEMAddElement(dm_contact, simple->getDomainFEName().c_str());
CHKERR DMMoFEMAddElement(dm_contact, simple->getBoundaryFEName().c_str());
CHKERR DMMoFEMSetSquareProblem(dm_contact, PETSC_TRUE);
CHKERR DMMoFEMSetIsPartitioned(dm_contact, PETSC_TRUE);
CHKERR DMSetUp(dm_contact);
MoFEMFunctionReturn(0);
}
in the main, this is how I set them:
auto dm = simple->getDM();
if (data.with_contact && data.with_plasticity) {
CHKERR set_contact_dm(mField, dm, dmContact);
CHKERR set_plastic_ep_dm(mField, dmContact, dmEP);
CHKERR set_plastic_tau_dm(mField, dmEP, dmTAU);
} else if (data.with_plasticity) {
CHKERR set_plastic_ep_dm(mField, dm, dmEP);
CHKERR set_plastic_tau_dm(mField, dmEP, dmTAU);
}
Lukasz Kaczmraczyk
Dec 7, 2021, 11:10:01 AM12/7/21
to MoFEM Q&A
Creation of sub DMs looks OK. What you mean that issue is for assembly to "EP". You have written this:
"When I apply a single FS contact (additive or multiplicative), it works. The second level with multiplicative (EP) - it works as well,
Second level with Schur - problems with assembly."
It stops to work when you add third level for TAU, EP stops to work?
Lukasz
When I apply a single FS contact (additive or multiplicative), it works. The second level with multiplicative (EP) - it works as well,
Lukasz Kaczmraczyk
Dec 7, 2021, 11:15:44 AM12/7/21
to MoFEM Q&A
Kariol,
Do you using "is_map" on assembly level? Could be that you change it?
Also, do you know assembly of what field to what matrix causing problem?
Lukasz
Reply all
Reply to author
Forward
0 new messages