> On 21. Sep 2017, at 17:33, Matthew Knepley <[email protected]> wrote:
>
> On Thu, Sep 21, 2017 at 11:28 AM, Maximilian Hartig <[email protected]
> <mailto:[email protected]>> wrote:
>> On 20. Sep 2017, at 22:51, Matthew Knepley <[email protected]
>> <mailto:[email protected]>> wrote:
>>
>> On Wed, Sep 20, 2017 at 3:46 PM, Maximilian Hartig <[email protected]
>> <mailto:[email protected]>> wrote:
>>
>>> On 20. Sep 2017, at 19:05, Matthew Knepley <[email protected]
>>> <mailto:[email protected]>> wrote:
>>>
>>> On Wed, Sep 20, 2017 at 12:57 PM, Maximilian Hartig
>>> <[email protected] <mailto:[email protected]>> wrote:
>>>> On 20. Sep 2017, at 18:17, Matthew Knepley <[email protected]
>>>> <mailto:[email protected]>> wrote:
>>>>
>>>> On Wed, Sep 20, 2017 at 11:46 AM, Maximilian Hartig
>>>> <[email protected] <mailto:[email protected]>> wrote:
>>>> Hello,
>>>>
>>>> I’m trying to implement plasticity using petscFE but I am quite stuck
>>>> since a while. Here’s what I’m trying to do:
>>>>
>>>> I have a TS which solves the following equation:
>>>> gradient(stress) +Forces = density*acceleration
>>>> where at the moment stress is a linear function of the strain and hence
>>>> the gradient of the displacement. This works fine. Now I want to compare
>>>> the stress to a reference value and if it lies above this yield stress, I
>>>> have to reevaluate the stress at the respective location. Then I need to
>>>> update the plastic strain / yield stress at this location.
>>>> I tried doing that first by solving three fields at the same time:
>>>> displacements, stresses and yield stress. This failed.
>>>> Then, I tried solving only for displacement increments, storing the
>>>> displacements, stresses and yield stress from the past time step in an
>>>> auxiliary field. The auxiliary fields are updated after each time step
>>>> with a second SNES, using the displacement increments from the current,
>>>> converged time step. This also failed.
>>>> In both cases the code had problems converging and when it did, I ended up
>>>> with negative plastic strain. This is not possible and I don’t know how it
>>>> happens because I explicitly only increment the plastic strain when the
>>>> increment is positive.
>>>>
>>>> I am sure there is an easy solution to how I can update the internal
>>>> variables and determine the correct stress for the residual but I just
>>>> cannot figure it out. I’d be thankful for any hints.
>>>>
>>>> It looks like there are two problems above:
>>>>
>>>> 1) Convergence
>>>>
>>>> For any convergence question, we at minimum need to see the output of
>>>>
>>>> -snes_view -snes_converged_reason -snes_monitor
>>>> -ksp_monitor_true_residual -snes_linesearch_monitor
>>>>
>>>> However, this does not seem to be the main issue.
>>>>
>>>> 2) Negative plastic strain
>>>
>>> This is what I’m mainly concerned with.
>>>>
>>>> If the system really converged (I cannot tell without other information),
>>>> then the system formulation is wrong. Of course, its
>>>> really easy to check by just plugging your solution into the residual
>>>> function too. I do not understand your explanation above
>>>> completely however. Do you solve for the plastic strain or the increment?
>>>
>>> I am trying to find a formulation that works and I think there is a core
>>> concept I am just not “getting”.
>>> I want to solve for the displacements.
>>> This works fine in an elastic case. When plasticity is involved, I need to
>>> determine the actual stress for my residual evaluation and I have not found
>>> a way to do that.
>>> All formulations for stress I found in literature use strain increments so
>>> I tried to just solve for increments each timestep and then add them
>>> together in tspoststep. But I still need to somehow evaluate the stress for
>>> my displacement increment residuals. So currently, I have auxiliary fields
>>> with the stress and the plastic strain.
>>>
>>> First question: Don't you get stress by just applying a local operator,
>>> rather than a solve?
>> That depends on the type of plasticity.
>>
>> What type of plasticity is not local?
> I did not express myself correctly. The plasticity is local, but for
> nonlinear hardening I would have to iteratively solve to get the correct
> stress and plastic strain from the displacement increments.
>
> Okay, I see. Lets leave that for the time being.
>>
>> For a linear hardening formulation it is correct that I could just apply a
>> local operator. I’d be happy with that for now. But I’d still need to save
>> stress state and plastic strain to determine whether or not I’m dealing with
>> a plasticity. I don’t know how to do that inside the residual evaluation.
>>
>> I do not know what you mean by this, meaning why you can't just save these
>> as auxiliary fields. Also, it would seem to be enough to have the old
>> displacement and the plastic strain.
> Yes, I can update the auxiliary fields but only after I solved for the
> displacements in my understanding. I still need the stress though as I have
> to determine whether I am in the plastic or the elastic domain.
>
> Right, but if the stress is local, then you can just compute it in the
> kernel (this is what I do). Even if the relationship is implicit, you can
> compute that solve local to the (Gauss) point.
Just to clarify, with kernel you mean the residual function I give to
PetscDSSetResidual(…), correct? I currently do this but did not find to save
plastic strain and stress inside this function. So I decided to reevaluate them
after I had the correct displacement field. Since I’d have to do that only once
per time step, I thought the extra computational effort would be negligible if
I could just project them.
>> Plus DMProjectField seems to have problems evaluating the gradient when
>> boundary conditions are imposed.
>>
>> There are several examples where we do exactly this. Can you show me what
>> you mean by this?
>
> Yes, of course. I sent an example a week ago but maybe there was a problem
> with the attached files. I’ll copy and paste the code and the gmsh file as
> text below.
>
> Okay, I got this. Can you tell me how to run it and why you think stuff is
> wrong?
Sure, I compile with the following makefile (evidently the c file is called
projectfieldtest.c):
PETSC_DIR = ...
PETSC_ARCH = arch-darwin-c-debug
CFLAGS =
FFLAGS =
CPPFLAGS =
FPPFLAGS =
MANSEC =
LOCDIR =
include ${PETSC_DIR}/${PETSC_ARCH}/lib/petsc/conf/petscvariables
include ${PETSC_DIR}/lib/petsc/conf/variables
include ${PETSC_DIR}/lib/petsc/conf/rules
projectfieldtest: projectfieldtest.o chkopts
-${CLINKER} -o projectfieldtest projectfieldtest.o ${PETSC_LIB}
${RM} projectfieldtest.o
Then run with the following options:
-disp_petscspace_order 2 -stress_petscspace_order 2
I use petsc from bitbucket origin/master with configure options:
--download-mumps --download-chacoa --download-scalapack --download-mpich
--with-fc=gfortran-mp-6 --download-triangle --with-cc=gcc-mp-6
--with-cxx=g++-mp-6 --download-fblaspack --download-parmetis --download-metis
I think the result is wrong because I compare it to a solution that I
calculated by hand. I project a displacement function to the dm, then compare
the stress I get from DMProjectField to the stress I calculated by hand and
then projected using DMProjectFunction.
Without this line:
ierr = DMAddBoundary(dm, PETSC_TRUE, "fixed", "Face
Sets",0,Ncomp,restrictall,(void (*)()) zero_vector,
Nfid,fid,NULL);CHKERRQ(ierr);
the expected solution is the same as the projected solution. With it, they
differ.
Thanks,
Max
>
> Thanks,
>
> Matt
>
> Thanks,
> Max
>
>
>
>
> #include <petscdm.h>
> #include <petscdmlabel.h>
> #include <petscds.h>
> #include <petscdmplex.h>
> #include <petscksp.h>
> #include <petscsnes.h>
> #include <petscts.h>
> #include <math.h>
> #include <petscsys.h>
>
> /* define the pointwise functions we'd like to project */
>
> void projectstress(PetscInt dim, PetscInt Nf, PetscInt NfAux, const PetscInt
> uOff[],
> const PetscInt uOff_x[], const PetscScalar u[],
> const PetscScalar u_t[],
> const PetscScalar u_x[], const PetscInt aOff[],
> const PetscInt aOff_x[],
> const PetscScalar a[], const PetscScalar a_t[],
> const PetscScalar a_x[],
> PetscReal t, const PetscScalar x[], PetscInt
> numConstants,
> const PetscScalar constants[], PetscScalar v[]){
> const PetscReal mu =76.923076923, lbda=115.384615385;
> PetscInt Ncomp = dim;
> PetscInt comp,d;
> PetscReal sigma[dim*dim];
>
> for(comp=0;comp<Ncomp;comp++){
> for(d=0;d<dim;d++){
> sigma[comp*dim+d]=mu*(u_x[comp*dim+d]+u_x[d*dim+comp]);
> }
> for(d=0;d<dim;d++){
> sigma[comp*dim+comp]+=lbda*u_x[d*dim+d];
> }
> }
>
> for(d=0;d<dim;d++){
> v[d] = sigma[d*dim+d];
> }
> v[3] = sigma[0*dim+1];
> v[4] = sigma[0*dim+2];
> v[5] = sigma[1*dim+2];
>
>
> }
>
> void projectdisplacement(PetscInt dim, PetscInt Nf, PetscInt NfAux, const
> PetscInt uOff[],
> const PetscInt uOff_x[], const PetscScalar u[],
> const PetscScalar u_t[],
> const PetscScalar u_x[], const PetscInt aOff[],
> const PetscInt aOff_x[],
> const PetscScalar a[], const PetscScalar a_t[],
> const PetscScalar a_x[],
> PetscReal t, const PetscScalar x[], PetscInt
> numConstants,
> const PetscScalar constants[], PetscScalar v[]){
> PetscInt d;
>
> for(d=0;d<dim;d++){
> v[d] =u[d];
> }
> }
>
> static PetscErrorCode initial_displacement_vector(PetscInt dim, PetscReal
> time, const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
> {
> u[0]=0.0;
> u[1]=0.1*pow(x[2],2);
> u[2]=0.1*x[2];
> return 0;
> }
>
> static PetscErrorCode expected_stress_vector(PetscInt dim, PetscReal time,
> const PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
> {
> const PetscReal mu =76.923076923, lbda=115.384615385;
> PetscReal strain[dim*dim];
> PetscReal gradu[dim*dim];
> PetscReal stress[dim*dim];
> PetscInt i,j;
>
> /* gradient of the displacement field: */
> for(i=0;i<dim;i++){
> for(j=0;j<dim;j++){
> gradu[i*dim+j]=0.0;
> }
> }
>
> gradu[1*dim+2]=0.2*x[2];
> gradu[2*dim+2]=0.1;
>
> for(i=0;i<dim;i++){
> for(j=0;j<dim;j++){
> strain[i*dim+j]=0.5*(gradu[i*dim+j]+gradu[j*dim+i]);
> }
> }
>
> for(i=0;i<dim;i++){
> for(j=0;j<dim;j++){
> stress[i*dim+j]=2.0*mu*strain[i*dim+j];
> }
> for(j=0;j<dim;j++){
> stress[i*dim+i]+=lbda*strain[j*dim+j];
> }
> }
>
> for(i=0;i<dim;i++){
> u[i] = stress[i*dim+i];
> }
> u[3] = stress[0*dim+1];
> u[4] = stress[0*dim+2];
> u[5] = stress[1*dim+2];
>
>
> return 0;
> }
>
> static PetscErrorCode zero_stress(PetscInt dim, PetscReal time, const
> PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
> {
> const PetscInt Ncomp = 2*dim;
> PetscInt comp;
> for (comp=0;comp<Ncomp;comp++) u[comp]=0.0;
> return 0;
> }
>
> static PetscErrorCode zero_vector(PetscInt dim, PetscReal time, const
> PetscReal x[], PetscInt Nf, PetscScalar *u, void *ctx)
> {
> PetscInt comp;
> for (comp=0;comp<dim;comp++) u[comp]=0.0;
> return 0;
> }
>
> static PetscErrorCode SetupDiscretization(DM dm){
>
> PetscInt dim = 3;
> PetscFE fe_displacement, fe_stress;
> PetscDS prob;
> PetscErrorCode ierr;
> PetscBool simplex = PETSC_TRUE;
> PetscQuadrature q;
> PetscInt order;
>
> /* get the dimension of the problem from the DM */
>
> ierr = DMGetDimension(dm, &dim); CHKERRQ(ierr);
>
> /* Creating the FE for the displacement */
> ierr = PetscFECreateDefault(dm, dim, dim,
> simplex,"disp_",PETSC_DEFAULT,&fe_displacement);CHKERRQ(ierr);
> ierr = PetscObjectSetName((PetscObject) fe_displacement,
> "displacement");CHKERRQ(ierr);
> ierr = PetscFEGetQuadrature(fe_displacement,&q);CHKERRQ(ierr);
> ierr = PetscQuadratureGetOrder(q,&order);CHKERRQ(ierr);
> /* Creating the FE for the stress */
> ierr =
> PetscFECreateDefault(dm,dim,2*dim,simplex,"stress_",PETSC_DEFAULT,&fe_stress);CHKERRQ(ierr);
> ierr = PetscFESetQuadrature(fe_stress,q);CHKERRQ(ierr);
> ierr = PetscObjectSetName((PetscObject) fe_stress,
> "cauchy_stress");CHKERRQ(ierr);
>
>
> /* Discretization and boundary conditons: */
> ierr = DMGetDS(dm, &prob);CHKERRQ(ierr);
> ierr = PetscDSSetDiscretization(prob, 0, (PetscObject) fe_displacement);
> CHKERRQ(ierr);
> ierr = PetscDSSetDiscretization(prob, 1, (PetscObject)
> fe_stress);CHKERRQ(ierr);
> ierr = DMSetDS(dm, prob); CHKERRQ(ierr);
>
> /* Define the boundaries */
>
> const PetscInt Ncomp = dim;
> const PetscInt Nfid = 1;
> PetscInt fid[Nfid]; /* fixed faces [numer of fixed faces] */
>
> PetscInt restrictall[3] = {0, 1, 2}; /* restricting all movement */
>
>
>
> fid[0] = 3; /* fixed face */
>
> ierr = DMAddBoundary(dm, PETSC_TRUE, "fixed", "Face
> Sets",0,Ncomp,restrictall,(void (*)()) zero_vector,
> Nfid,fid,NULL);CHKERRQ(ierr);
>
>
> ierr = PetscFEDestroy(&fe_displacement); CHKERRQ(ierr);
>
> ierr = PetscFEDestroy(&fe_stress); CHKERRQ(ierr);
>
> return(0);
> }
> int main(int argc, char *argv[])
> {
> DM dm, distributeddm; /* problem definition */
> Vec u,expected_solution,projected_solution;
> PetscViewer viewer;
> int dim; /* dimension of the anlysis */
> PetscErrorCode ierr;
> PetscMPIInt rank, numProcs;
> PetscPartitioner part;
>
>
>
> ierr = PetscInitialize(&argc,&argv,NULL,NULL);
> ierr = DMPlexCreateFromFile(PETSC_COMM_WORLD,"testcube.msh",
> PETSC_TRUE,&(dm));
> ierr = DMGetDimension(dm,&(dim));
>
>
> /* distribute the mesh */
> MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
> MPI_Comm_size(PETSC_COMM_WORLD, &numProcs);
> DMPlexGetPartitioner(dm, &part);
> PetscPartitionerSetType(part, PETSCPARTITIONERPARMETIS);
>
> ierr = DMPlexDistribute(dm,0,NULL,&distributeddm);
>
> if (distributeddm) {
> ierr=DMDestroy(&(dm));
> dm = distributeddm;
> }
>
>
> ierr = DMView(dm,PETSC_VIEWER_STDOUT_WORLD);
> ierr = DMSetMatType(dm, MATAIJ);CHKERRQ(ierr);
>
>
>
> ierr = SetupDiscretization(dm);CHKERRQ(ierr);
>
> ierr = DMPlexCreateClosureIndex(dm,NULL);CHKERRQ(ierr);
> ierr = DMCreateGlobalVector(dm, &(u)); CHKERRQ(ierr);
> ierr = VecDuplicate(u,&expected_solution); CHKERRQ(ierr);
> ierr = VecDuplicate(u,&projected_solution); CHKERRQ(ierr);
>
> /* intitialize the fields: */
>
> PetscErrorCode (*initial[2])(PetscInt dim, PetscReal time, const PetscReal
> x[], PetscInt Nf, PetscScalar u[], void* ctx) =
> {initial_displacement_vector,zero_stress};
> ierr = DMProjectFunction(dm,0.0,initial,NULL, INSERT_ALL_VALUES, u);
> CHKERRQ(ierr);
>
> PetscErrorCode (*expected_sol[2])(PetscInt dim, PetscReal time, const
> PetscReal x[], PetscInt Nf, PetscScalar u[], void* ctx) =
> {initial_displacement_vector,expected_stress_vector};
> ierr = DMProjectFunction(dm,0.0,expected_sol,NULL, INSERT_ALL_VALUES,
> expected_solution); CHKERRQ(ierr);
>
> ierr =
> PetscViewerVTKOpen(PETSC_COMM_WORLD,"expected_solution.vtu",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
> ierr = PetscObjectSetName((PetscObject) expected_solution,
> "expected_fields"); CHKERRQ(ierr);
> ierr = VecView(expected_solution,viewer);CHKERRQ(ierr);
> ierr= PetscViewerDestroy(&viewer); CHKERRQ(ierr);
>
>
>
> /* project the fields: */
>
> void (*projection[2])(PetscInt dim, PetscInt Nf, PetscInt NfAux, const
> PetscInt uOff[],
> const PetscInt uOff_x[], const PetscScalar u[],
> const PetscScalar u_t[],
> const PetscScalar u_x[], const PetscInt aOff[],
> const PetscInt aOff_x[],
> const PetscScalar a[], const PetscScalar a_t[],
> const PetscScalar a_x[],
> PetscReal t, const PetscScalar x[], PetscInt
> numConstants,
> const PetscScalar constants[], PetscScalar v[]) =
> {projectdisplacement, projectstress};
>
>
> ierr = DMProjectField(dm, 0.0, u,
> projection,INSERT_ALL_VALUES,projected_solution); CHKERRQ(ierr);
> ierr =
> PetscViewerVTKOpen(PETSC_COMM_WORLD,"projected_solution.vtu",FILE_MODE_WRITE,&viewer);CHKERRQ(ierr);
> ierr = PetscObjectSetName((PetscObject) projected_solution,
> "projected_fields"); CHKERRQ(ierr);
> ierr = VecView(projected_solution,viewer);CHKERRQ(ierr);
> ierr= PetscViewerDestroy(&viewer); CHKERRQ(ierr);
>
>
> VecDestroy(&u);
> VecDestroy(&expected_solution);
> VecDestroy(&projected_solution);
> DMDestroy(&dm);
>
> PetscFinalize();
>
> return 0;
> }
>
>
> testcube.msh:
>
> $MeshFormat
> 2.2 0 8
> $EndMeshFormat
> $PhysicalNames
> 6
> 2 1 "back"
> 2 2 "front"
> 2 3 "bottom"
> 2 4 "right"
> 2 5 "top"
> 2 6 "left"
> $EndPhysicalNames
> $Nodes
> 10
> 1 0 -0.05 0
> 2 0.1 -0.05 0
> 3 0.1 0.05 0
> 4 0 0.05 0
> 5 0 -0.05 0.1
> 6 0.1 -0.05 0.1
> 7 0.1 0.05 0.1
> 8 0 0.05 0.1
> 9 0.05 0 0
> 10 0.05 0 0.1
> $EndNodes
> $Elements
> 28
> 1 2 2 1 6 1 2 9
> 2 2 2 1 6 1 9 4
> 3 2 2 1 6 2 3 9
> 4 2 2 1 6 3 4 9
> 5 2 2 3 15 5 1 6
> 6 2 2 3 15 1 2 6
> 7 2 2 4 19 6 2 7
> 8 2 2 4 19 2 3 7
> 9 2 2 5 23 7 3 8
> 10 2 2 5 23 3 4 8
> 11 2 2 6 27 8 4 1
> 12 2 2 6 27 8 1 5
> 13 2 2 2 28 5 6 10
> 14 2 2 2 28 5 10 8
> 15 2 2 2 28 6 7 10
> 16 2 2 2 28 7 8 10
> 17 4 2 29 1 1 2 9 6
> 18 4 2 29 1 6 5 10 9
> 19 4 2 29 1 1 5 6 9
> 20 4 2 29 1 1 9 4 8
> 21 4 2 29 1 10 5 8 9
> 22 4 2 29 1 9 5 8 1
> 23 4 2 29 1 2 3 9 7
> 24 4 2 29 1 7 6 10 9
> 25 4 2 29 1 2 6 7 9
> 26 4 2 29 1 3 4 9 8
> 27 4 2 29 1 8 7 10 9
> 28 4 2 29 1 3 7 8 9
> $EndElements
>
>
>
>
>>
>> Thanks,
>>
>> Matt
>>
>> Thanks,
>> Max
>>>
>>> Thanks,
>>>
>>> Matt
>>>
>>> I evaluate the current trial stress by adding a stress increment assuming
>>> elastic behaviour. If the trial stress lies beyond the yield stress I
>>> calculate the corrected stress to evaluate my residual for the
>>> displacements. But now I somehow need to update my plastic strain and the
>>> stress in the auxiliary fields. So in tspoststep I created another SNES to
>>> now calculate the stress and plastic strain while the displacement is the
>>> auxiliary field.
>>>
>>> I’m sure there’s an elegant solution on how to update internal variables
>>> but I have not found it.
>>>
>>> Thanks,
>>> Max
>>>>
>>>> Thanks,
>>>>
>>>> Matt
>>>>
>>>> Thanks,
>>>> Max
>>>>
>>>>
>>>>
>>>> --
>>>> What most experimenters take for granted before they begin their
>>>> experiments is infinitely more interesting than any results to which their
>>>> experiments lead.
>>>> -- Norbert Wiener
>>>>
>>>> http://www.caam.rice.edu/~mk51/
>>>
>>>
>>>
>>>
>>> --
>>> What most experimenters take for granted before they begin their
>>> experiments is infinitely more interesting than any results to which their
>>> experiments lead.
>>> -- Norbert Wiener
>>>
>>> http://www.caam.rice.edu/~mk51/
>>
>>
>>
>>
>> --
>> What most experimenters take for granted before they begin their experiments
>> is infinitely more interesting than any results to which their experiments
>> lead.
>> -- Norbert Wiener
>>
>> http://www.caam.rice.edu/~mk51/
>
>
>
>
> --
> What most experimenters take for granted before they begin their experiments
> is infinitely more interesting than any results to which their experiments
> lead.
> -- Norbert Wiener
>
> http://www.caam.rice.edu/~mk51/
