> On 20. Sep 2017, at 22:51, Matthew Knepley <knep...@gmail.com> wrote: > > On Wed, Sep 20, 2017 at 3:46 PM, Maximilian Hartig <imilian.har...@gmail.com > <mailto:imilian.har...@gmail.com>> wrote: > >> On 20. Sep 2017, at 19:05, Matthew Knepley <knep...@gmail.com >> <mailto:knep...@gmail.com>> wrote: >> >> On Wed, Sep 20, 2017 at 12:57 PM, Maximilian Hartig >> <imilian.har...@gmail.com <mailto:imilian.har...@gmail.com>> wrote: >>> On 20. Sep 2017, at 18:17, Matthew Knepley <knep...@gmail.com >>> <mailto:knep...@gmail.com>> wrote: >>> >>> On Wed, Sep 20, 2017 at 11:46 AM, Maximilian Hartig >>> <imilian.har...@gmail.com <mailto:imilian.har...@gmail.com>> 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. > > 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. > > 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. 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/ <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/ <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/ <http://www.caam.rice.edu/~mk51/>