Sorry, please disregard my last question. I've reviewed the Step-20 solver that I tried to use and I see that it is a different set-up to my problem. I'm now working to apply the Schur complement in the correct manner.
On Wednesday, November 8, 2023 at 10:32:35 AM UTC-5 Alex Quinlan wrote: > Alright, I've had some decent success, though I've hit another snag in the > solver section. Your suggestions have opened the doors and actually > allowed me to shift gears to my preferred approach. > > I've changed the end of Setup_System() to be: > > dsp.collect_sizes(); > > *DoFTools::make_sparsity_pattern(dof_handler_sld, dsp.block(0,0)); > DoFTools::make_sparsity_pattern(dof_handler_sh, dsp.block(1,1));* > > sparsity_pattern.copy_from(dsp); > system_matrix.reinit(sparsity_pattern); > system_rhs.reinit(dofs_per_block); > solution.reinit(dofs_per_block); > std::cout << " " << std::endl; > > Doing so allowed me to successfully run your suggested code: > > for (unsigned int i = 0; i < dofs_per_cell_sld; ++i) > for (unsigned int j = 0; j < dofs_per_cell_sld; ++j) > system_matrix.*block(0,0)*.add(local_dof_indices_sld[i], > > local_dof_indices_sld[j], > cell_matrix_sld(i, j)); > > > This success prompted me to try a different approach to applying the cell > stiffness matrix to the global stiffness matrix. I have an object of > AffineConstraints, so I wanted to use the > AffineConstraints::distribute_local_to_global(). > > For the block system, I prepared the constraints like this: (Note: > nodal_bcs() is a routine that I've written to apply constraints that I've > read in from a .csv file) > > // Set up Solid Constraints > constraints_sld.clear(); > nodal_bcs<3,3>(solidbcfile, dof_handler_sld ,constraints_sld); > constraints_sld.close(); > > // Set-up Shell Constraints > constraints_sh.clear(); > nodal_bcs<2,3>(shellbcfile, dof_handler_sh ,constraints_sh); > constraints_sh.close(); > > // Combine Solid and Shell constraints > constraints.copy_from(constraints_sld); > constraints_sh.shift(n_dof_sld); > constraints.merge(constraints_sh ); > > This allowed me to add this line to the Cell loop in the Assemble_System(): > > constraints.distribute_local_to_global( > cell_matrix_sld, cell_rhs_sld, local_dof_indices_sld, > system_matrix.block(0,0), system_rhs.block(0)); > > > All of this now seems to work and the matrices seem to be correct. > > ----- > > Now, I'm having trouble with the solver. I think the cause is the empty > off-diagonal matrix I am using the approach outlined in Step-20: > > const auto &M = system_matrix.block(0, 0); > const auto &B = system_matrix.block(0, 1); > const auto &F = system_rhs.block(0); > const auto &G = system_rhs.block(1); > auto &U = solution.block(0); > auto &P = solution.block(1); > const auto op_M = linear_operator(M); > const auto op_B = linear_operator(B); > > ReductionControl reduction_control_M(2000, 1.0e-18, 1.0e-10); > SolverCG<Vector<double>> solver_M(reduction_control_M); > PreconditionJacobi<SparseMatrix<double>> preconditioner_M; > preconditioner_M.initialize(M); > > const auto op_M_inv = inverse_operator(op_M, solver_M, > preconditioner_M); > const auto op_S = transpose_operator(op_B) * op_M_inv * op_B; > const auto op_aS = transpose_operator(op_B) * > linear_operator(preconditioner_M) * op_B; > > IterationNumberControl iteration_number_control_aS(30, 1.e-18); > SolverCG<Vector<double>> solver_aS(iteration_number_control_aS); > const auto preconditioner_S = inverse_operator(op_aS, solver_aS, > PreconditionIdentity()); > const auto schur_rhs = transpose_operator(op_B) * op_M_inv * F - G; > SolverControl solver_control_S(2000, 1.e-12); > SolverCG<Vector<double>> solver_S(solver_control_S); > const auto op_S_inv = inverse_operator(op_S, solver_S, > preconditioner_S); > > * P = op_S_inv * schur_rhs;* > std::cout << solver_control_S.last_step() > << " CG Schur complement iterations to obtain convergence." > << std::endl; > U = op_M_inv * (F - op_B * P); > > At the red line above, I am given this error: > > An error occurred in line <637> of file </home/.../lac/solver_cg.h> in > function > void dealii::internal::SolverCG::IterationWorker<VectorType, > MatrixType, PreconditionerType, <template-parameter-1-4> > >::do_iteration(unsigned int) [with VectorType = dealii::Vector<double>; > MatrixType = dealii::LinearOperator<dealii::Vector<double>, > dealii::Vector<double>, > dealii::internal::LinearOperatorImplementation::EmptyPayload>; > PreconditionerType = dealii::LinearOperator<dealii::Vector<double>, > dealii::Vector<double>, > dealii::internal::LinearOperatorImplementation::EmptyPayload>; > <template-parameter-1-4> = int] > The violated condition was: > std::abs(p_dot_A_dot_p) != 0. > Additional information: > A piece of code is attempting a division by zero. This is likely going > to lead to results that make no sense. > > I'm guessing this is related to the inverse operators? So, I'm curious if > there's some way to use this approach for problems that have a Schur > complements that can be either zero or non-zero. > > Thanks, > Alex > > On Wednesday, November 8, 2023 at 9:50:17 AM UTC-5 Alex Quinlan wrote: > >> Dear Luca and Wolfgang, >> >> Thanks for your pointers; I will work to implement them. >> >> My github account is: https://github.com/AlexQuinlan >> >> >> Thanks, >> Alex >> >> On Wednesday, November 8, 2023 at 3:36:43 AM UTC-5 luca....@gmail.com >> wrote: >> >>> Dear Alex, >>> >>> if you send me your github user, I’ll give you access to a code that >>> does exactly what you are trying to do, so that you can get some >>> inspiration on how to assemble the coupling matrices. >>> >>> There is an open PR (the review has been stalling for some time on that) >>> https://github.com/dealii/dealii/pull/15773 that tackle exactly this >>> use case. >>> >>> Let me know if you need some assistance. >>> >>> L. >>> >>> > On 8 Nov 2023, at 05:02, Wolfgang Bangerth <bang...@colostate.edu> >>> wrote: >>> > >>> > >>> > Alex: >>> > >>> > > [...] >>> > >>> > I think everything up to this point looks reasonable. >>> > >>> >> From what I can see, the system_matrix looks right with n_blocks = 2 >>> and start_indices = {0, 24, 72}, as expected. >>> >> Next, I run through the Assemble_System() to create a local cell >>> matrix. I start with the Solid elements and cell_matrix_sld. I'm skipping >>> the calculation of the cell matrix for brevity's sake. >>> >> for (const auto &cell : dof_handler_sld.active_cell_iterators()) >>> >> { >>> >> // *** Placeholder for calculating cell_matrix_sld, a 24x24 matrix >>> *** >>> >> cell->get_dof_indices(local_dof_indices_sld); >>> >> // Add the cell matrix to the system matrix >>> >> for (unsigned int i = 0; i < dofs_per_cell_sld; ++i) >>> >> for (unsigned int j = 0; j < dofs_per_cell_sld; ++j) >>> >> system_matrix.add(local_dof_indices_sld[i], >>> >> local_dof_indices_sld[j], >>> >> cell_matrix_sld(i, j)); >>> >> } // end of cell loop >>> >> The bit of code 'system_matrix.add(...)' works for when i=0 and j=0, >>> but fails on the next loop when j=1. I get this error: >>> > >>> > The issue here is that you're indexing into the *global* matrix with >>> the indices of only one of the DoFHandler objects. That can't go well. It >>> *should* work if you do the indexing only in the matrix block that >>> corresponds to dof_handler_sld: >>> > >>> > system_matrix.block(block_index_sld, block_index_sld) >>> > .add(local_dof_indices_sld[i], >>> > local_dof_indices_sld[j], >>> > cell_matrix_sld(i, j)); >>> > >>> > where you define block_index_sld somewhere to either be zero or one. >>> > >>> > In your case, this won't work yet because you hadn't built the >>> sparsity pattern yet. You can do that the same way: >>> > >>> > DoFTools::make_sparsity_pattern(dof_handler_sld, >>> > block_dsp.block(block_index_sld, block_index_sld)); >>> > >>> > Best >>> > W. >>> > >>> > -- >>> > >>> ------------------------------------------------------------------------ >>> > Wolfgang Bangerth email: bang...@colostate.edu >>> > www: http://www.math.colostate.edu/~bangerth/ >>> > >>> > >>> > -- >>> > The deal.II project is located at http://www.dealii.org/ >>> > For mailing list/forum options, see >>> https://groups.google.com/d/forum/dealii?hl=en >>> > --- You received this message because you are subscribed to the Google >>> Groups "deal.II User Group" group. >>> > To unsubscribe from this group and stop receiving emails from it, send >>> an email to dealii+un...@googlegroups.com. >>> > To view this discussion on the web visit >>> https://groups.google.com/d/msgid/dealii/9a890d55-91f6-ec16-98be-acdb730919d3%40colostate.edu. >>> >>> >>> >>> -- The deal.II project is located at http://www.dealii.org/ For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en --- You received this message because you are subscribed to the Google Groups "deal.II User Group" group. To unsubscribe from this group and stop receiving emails from it, send an email to dealii+unsubscr...@googlegroups.com. To view this discussion on the web visit https://groups.google.com/d/msgid/dealii/50795d91-1f5f-42e2-85f5-4f6e7fca353cn%40googlegroups.com.
